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ABSTRACT 


The filtering problem is studied for a system compo- 
sed of an inertial navigation system giving continuous in- 
dication of position and velocity, and a radar or some 
other external device giving continuous or discrete posi- 
tional information. After a brief review of the results of 
the filtering theory using the state variables representa- 
tion, an error analysis yields three different mathematical 
models for the I.N.S., represented by a set of linear diffe- 
rential equations that can be written with the state vari- 
ables method. Prom there, the equations for both the conti- 
nuous and the discrete filtering schemes are derived, assu- 
ming the only error sources are a white noise at the acce- 
leration level in the I.N.S. and a white or Gaussian (in 
the discrete measurement mode) noise in the radar, functio- 
nal relations are obtained between the position and veloci- 
ty root mean square errors and the following characteristic 
quantities : noises in both the I.N.S. and the radar and 0 - 
perating time (time between the measurements) for the dis- 
drete filter scheme. 
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NOTATIONS 


Underlined lower or upper case letters are column 
vectors. For instance x denotes a vector with, components: 

x^, x 2 , ,x n ; x' is the transpose of x and is a row 

vector. Capital letters denote matrices, except for the 
noises N and R which can he either scalar or matrices. The 
transpose of the matrix M is denoted M' . 


Scalars 

t time 

t“ time just before the measurement at t„ 

time just after the measurement at t ffi 
m 

R earth radius 

e 

^ earth rate 
le 

C s , Cy,C z misalignment angles of the platform with 
respect to the navigation frame 

1 latitude 

1 longitude 

th. 

p. . elements of the matrix P ; i row, j column 

x 

m. . continuous filtering compensation terms 
a j 

th th 

k- - elements of the gain matrix K; i column, j row 
^ j 

a, h, c parameters for models 1 and 2 
RMS - RMY r.m.s. north and east position errors 

RVX - RVX r.m.s. north and east velocity errors 

OPT or At operating time 


Vectors 

x state- vector 

u I.N.S. noise vector 

v radar noise vector 
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£ output of the linear system 
m measurement vector 

x optimum estimate of x 

f specific force 

£ gravity vector 

C misalignment angles vector 

—ah an S alar of frame b with respect to frame a 

coordinatized in frame c 


Matrices 

F system parameters matrix 

G- input distribution matrix 

H output distribution matrix 

<|> (t,s) state transition matrix from time s to time t 
I unit matrix 

Q or N I.N.S. noise power spectral densities matrix 
R radar noise power spectral densities matrix or co- 
variance matrix 
[PF] performance function 
[P J velocity computation matrix (page 23 ) 

K gain matrix 


Others 


p Laplace operator 

x time derivative of the quantity x 

c subscript denoting a variable computed by the system 
n subscript or superscript denoting navigation frame 

cm " " " platform frame 

i " " 11 inertial frame 

e " " " earth frame 

sk superscript denoting a skew-symmetric matrix 
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CHAPTER I 
INTRODUCTION 


The design of an optimal filter to yield actual posi- 
tion and velocity, when given data from an Inertial Naviga- 
tion System and radar, depends upon the model chosen to re- 
present the I.N.S. 

The degree of sophistication of the filter must be 
matched with the accuracy of the available instruments and 
measurements, and since a short operating time for the reset 
of t Re I.N.S. lowers the accuracy requirements, this fact 
must also be considered in designing the filter. 

We can consider 3 models for the I.N.S.: 
model I : platform misalignment and cross-coupling between 
the axes are neglected; so, it is possible to study the two 
channels separately. 

model 2 : platform misalignment angles are introduced but 
the cross-coupling is still neglected. 

model 3 : finally, both misalignment angles and cross cou- 
pling are taken into account; this is theoretically the most 
accurate filter, but also the most complex one. 

It will be shown later that the steady state r.m.s. 
errors do not depend on the model chosen to represent the 


I.N.S 
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Although the filter cannot work in a continuous way 
because of the time required for the computations of the 
varying gains, it may appear useful to study this continuous 
filtering which is the optimum and may be used as a refe- 
rence for choosing the parameters of the discrete filter. 

As a first approach to this problem, only two error 
sources will be considered: I.N.S. noise which appears as 
a white noise at the acceleration level and measurement noi- 
se (radar noise) which is also considered as a white noise. 
These two noises are to be characterized in the following 
way: the I.R.S.noise by its power spectral density N or Q 
in all the' cases; the measurement noise by R, which is a 
power spectral density in the continuous process and a co- 
variance matrix in the discrete one. This difference will 
be necessary in part 5-3, in order to relate any discrete 
process to one particular continuous one. 

From both the theory and the results of the computa- 
tion, functional relations are derived between the noises, 
the operating time and the 'steady state' r.m.s. position 
and velocity errors. The final charts 17 through 19 enable 
one to choose the instruments in order to match the accura- 
cy requirements for a given mission. 

It is to be understood that , whenever the position ra- 
dar appears in the text, any other external position infor- 
mation can be used instead, including radio-navigation and 
observation, provided that the noise in this information 
is a Gaussian white noise. 
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Furthermore, only the errors in the optimum estimate 

, i 

are to be studied. This means that we shall not study the 
way this estimate must be generated from both informations. 
Although some equations as well as some signal flow dia- 
grams for this estimator are given in this paper, only the 
variance equation, yielding the covariance matrix and the 
errors, is solved. 
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CHAPTER 2 

FILTERING- THEORY : PRINCIPAL RESULTS 


None of the equations' of the filtering theory are de- 
rived. Only the principal results that will be used throug- 
hout this paper are summarized. Further information about 

these equations and their derivation can be found in refe- 

1 2 
rences and 

2,1 State transition method 

Any linear dynamical system described by a set of or- 
dinary differential equations can always be represented by 
the single equation: 

x(t) = F(t) x(t) + G(t) u(t) (2-1) 

x(t) is the state vector of dimension n , its coordi- 
nates x^ are the state variables. 

F( t ) is the system parameters matrix . 

G(t) is the (n.l) input distribution matrix . 

u(t) is a vector of dimension 1 called control vector . 

It is assumed that both matrices F(t) and G(t) are continu- 
ous functions of time. 
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The system is completely described when the vector 
output j£(t) f which can be of any dimension m less than n, 
is written as a linear combination of the state variables: 


£(t) « H(t ) x(t) (2-2) 

where H(t) is the (m*n) output distribution matrix . 

The block diagram of this system is given on figure 1 . 


P(t) represents the dynamics of the system. 

G(t) " the constraints due to inputs. 

H(t) " the constraints on observing the system 

from outputs. 

The general solution of this system (see ^ and is: 
x(t) = ^>(t,t 0 ) x(t Q ) +| G(s) u(s) da (2-3) 

K 

where <|>(t,s) is the state transition matrix , which is al- 
ways nonsingular and obeys the differential equation: 


c|>(t,s) 


= F(t) <|>(t,s) 


(2-4) 


with the "initial condition" (t,t) = I (unit matrix). 


2,2 Stochastic linear differential equation 

When a linear dynamical system is driven by some vec- 
tor noise v(t), its state obeys the equation (2-1), namely: 



(2-5) 
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x(t) = F(t) x(' fc ) + G(t) u(t) 


If u(t) is a white noise - i.e. a noise the power 
spectral density of which is a constant over some bandwidth- 
this equation becomes a stochastic linear differential e- 
quation and its solution is again 

G-(s) u( s ) ds 

In order to study the statistics of x(t), we define 


x(t) 


= <J> (t,t 0 ) x(t 0 ) + J <J)(t,s) 


mean of x(t) = m x (t) ~ E £x ( t )J 

covariance matrix of x(t) = P x (t) = E ^(t) .x 1 (t)J 
Then, it can be easily found that: 

m x (t) = <£(t,t e ) m x (t 0 ) 

P x (t) = $(t,t 0 -) P x (t 0 ) <£'(t,t 0 ) 

+ / <J>(t,s) G(s) Q(s) G'(s) <J>-'(t f s) ds 
% (2-6) 

where Q.(s) will be called the strength of the white noise 
and is defined by the relation 

E^u(t) .u' (t+s)] = Q(t) £ (s) (2-7) 

S (s) is a 1-dimension vector delta function such that 


£ (s) = 0 if s is not 0 

,+<© /+GO /-‘■CD 


,♦<30 / + 00 

/ ds-L / ds 2 / ds^ £ ( s ) — 1 

Leo '-<» L<o 
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2.3 The general filtering theory 

The message is a random process x(t) generated by a 
model obeying a stochastic linear . differential equation: 

x(t) = E(t) x(^) + G(t) u(t) (2-3) 

-The observed signal , or measurement is 

m(t) = H(t) x(t) + v(t) (2-9) 

where u(t) and v(t) are white noises with zero means and co- 
variance matrices 

? u (t) = E £u( t ) u'(s)j = Q(t) S(t-s) 

P v (t) = Ejj(t) v* (s)J = R(t) S(t-s) 
and Ejju(t) v(t)'J = 0 

We assume the matrices Q(t) and R(t) are non-negative defi- 
nite matrices continuously differentiable. 

figure 2 represents the block diagram of this system. 
The filtering problem is then to determine from the 
measurement m(t) the best estimate x(t), best in the sen- 
se that it maximizes the conditional probability density 
function fx/M^AO of the state vector x conditioned on 
the values (a priori past values as- well as actual ones) 
of the measurements. 

We only give the results of this theory: see ^ for 
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5 

free system and continuous filtering, and ' for discrete 
filtering. 


2.3.1 free system 

The system simply obeys the equation (2-8). No mea- 
surement is taken so that the best estimate £(t) is m^(t) 
with variance P (t) given by equation (2-6). 

A 

It is easier to write these equations as differential 

equations instead of integral ones. Differentiating (2-6) 

/ 

and using (2-4) yield: 


,x(t) = P(t) x(t) 

,P x (t) = P(t) P x (t) + P x (t) P'(t) + G(t) Q(t) G»(t) 

( 2 - 10 ) 


2.3.2 Continuous filtering 

Measurement m(t) is now taken in a continuous way. 
It can be shown that, in this case, the above equations 
become : 


x(t) = P(t)x(t) + P x (t)H:' (t)R 1 (t)[m(t)-H(t)x(t)J 
> = PP + PP' 4 - GQG* - PH , R’ 1 HP (2-11) 

A 

In the last equation the variable t has been dropped and 
this will stand throughout this paper whenever there is no 
ambiguity. 


2.3.3 Discrete filtering 

The measurement can only be taken at discrete times. 
The time between two measurements is assumed constant and 
called the operating time . The scheme is the following oner 
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Between the measurements, the system is free and obeys (2-10) 
At each measurement, we compute 

K( t ) = P»(t )H‘(t ) [k(t )P' (t )H* (t )+R(t )1 _1 
m' ' m m L m m ' m m J 

Then the best estimate is given by 

x(t m ) = K(t m )m(t m ) with covariance matrix: 

= [l-K(t n )H(t m )j P x (tp (2-12) 

In this last equation, ' means just before the measurement 
and ' + ' just after the measurement. 

The diagram of the continuous filter is given in 
figure 3 . 

2.4 The filtering problem in this paper 

The general filtering theory assumes that one can ta- 
ke measurements of some coordinates of the state vector x 
which is itself generated by the system. 

The problem is not quite the same here and is to be 
understood in a different way: 

1. the state vector x(t), which can include position and 
velocity informations obeys the homogeneous equation: 

x(t) = F(t) x(t) 

where P(t) depends primarily upon the trajectory. 

2. some indications about position, velocity and other co- 
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ordinates of the state vector are available as outputs of 
an inertial navigation system. But due to measurements and 
instruments inaccuracies, some noise is added so that the 
output obeys equation (2-8); we limit ourselves in this pa- 
per to the case where the only noise source is the gyro 
drift. Furthermore, following and' , we assume that this 
drift can be well enough approximated by a white noise at 
the acceleration level . 

3. on another hand, some components of the state vector are 
available as outputs to an external device - external to 
the inertial system - This includes altimeter, Doppler and 
position radar, etc... Here again, the inaccuracy yields 
a noise term in equation (2-9) which represents the exter- 
nal information. 

Thus, we have here two different measurements of the 
same vector and we want to get from them the best estimate 
of this vector. 

From now on, the external device will be a position 
radar. The measurement will be the difference between the 
position given by the inertial navigator and the position 
given by the radar. This measurement is considered as the' 
error in the indication of the position. Filtering this po- 
sition error yields the best estimate of the error on the 
state vector. The best estimate of the state vector itself 
is then obtained by substracting the estimated error from 
the output of the inertial system. 
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The general block diagram of these operations is gi- 
ven in figure 4 . 

This approach allows us to use an error analysis for 
the inertial navigator. . 

The following chapter investigates the different pos- 
sible models for this navigator. 
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CHAPTER 3 

THE THREE POSSIBLE MODELS 


3.1 Introduction 

THe complexity of- the filter increases with the com- 
plexity of the model representing the inertial .navigator. 
Indeed the filter gain K(t) depends on P(t) and the dimen- 
sions of the matrix E increase with the number of variables 
in the system. 

On another hand, the accuracy obtained after filtering 
depends both on the complexity of the model and on the "po- 
wer" of the noise. Thus it could be useless to filter with 
a complex model if the noise is important. 

In this paper, 3 models are analysed: 

1. the first one neglects both platform misalignment and 
cross-coupling between the axes. 

2. the second one takes into account the platform misalign- 
ment angles . . 

3. the last one is a complete 3-axes model- but assumes a 
steady state - i.e. "en route" conditions-. 

In this chapter, we determine for all of these models 
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the various matrices involved in the filtering theory: F, 
G, H, Q, R. 


3.2 Model 1 : no cross-coupling; no misalignment 


It is simply considered that the velocity vector is 
the time derivative of the position vector. A simple dia- 
gram of this I.N.S. is given in figure 5 . 

Consider the state vector x composed of 'the position 


and velocity errors as: 


x = 


dr. 

dr 


x 


y 


dv. 


x 


dv. 


y 

Then the equation (2-8) can be written as 


dr 

p 

.X 


dr 


*y 


dv/ 


.X 


dv 


y 

L 


x 


0 0 10 
0 0 0 1 
0 0 0 0 
0 0 0 0 

F 



'dr ] 


' 0 o' 


X 




dr 


0 0 


<i 

x{ «< 


1 0 


dv Tr 


0 1 

J 

L 




x 


G 


u 


x 


u. 


yj 


u 


(3-1-1) 


In this equation u x and u y are the two white noises on both 
the x and the y channels. 

The measurement equation (2-9) is: 


or 


dr 

X 


*i 

o' 


fdr 1 

X 


' v x 

dr 

«y 


0 

1 


dr 

y 

to ™ 

+ 

>. 

m 

— 

H 


X 

+ 

V 


(3-1-2) 
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Finally, the matrices Q and R are here: 



"n 

o’ 


H 

o’ 

Q = 

0 

N 

R = 

0 

R 


where N is the power spectral density of the I.N.S. noise 
and R " " 11 radar noise. 

3,3 Model 2 ; no cross-coupling 

The functional block diagram of a single axis iner- 
tial navigator instrumenting the navigation frame is given 
in figure 8 . The platform misalignment is now accounted 
for: this means that a component of the gravity vector is 
falsely sensed by at least one of the accelerometers. 

A misalignment angle about the z-axis (down) does 
not introduce a large error because the acceleration of the 
vehicle is usually small in comparison to g. On the contra- 
ry, a misalignment angle C (C ) about the x-axis (y-axis) 

x y 

causes the accelerometer along the y-axis (x-axis) to sen- 
se a component of gravity - g C„ (- g C„) in the small 

x y ' 

angles approximation . 

Thus we neglect the angle C about the down-axis; 

z 

then there is no coupling between the x( north) -axis and 
the y( east) -axis. The z( vertical) channel follows, exactly 
the equations of model 1 and the other two channels can be 
studied, separately. 

As shown in figure 9. the misalignment angles C__ and 
C are defined as correction angles -i.e. they are the an- 

*7 

gles the navigation frame should rotate about its X and Y 
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axes to align, itself with the platform axes 

Throughout this paper the accelerometer outputs are 
supposed to be the components of the specific force f along 
the sensitive axes of any of these accelerometers. The spe- 
cific force here is gravity minus acceleration ; so: 

£ = & - e‘i s S p 


The effects of C and C are the following: when C 

x y x 

is positive, the instrumented east axis is below the hori- 
zon; therefore the y-accelerometer senses + g G , When 

C is positive, the instrumented north axis is above the 

y 

horizon; therefore the x-accelerometer senses - g C . 

u 

Then, when the vehicle is moving, the accelerometers sense 


-R p L - g C along the x-axis 

e y * 

[-R cosl p a l + g C along the y-axis 


R q is the earth radius; L the latitude and 1 the longitude. 
© 


Therefore, the errors are 


S a 

0, on p L 

R y 
e 


S 

R 


C on cosh p 1 


The error model is given in figure 10 . 
We can write: 


P (dr x ) = dv x 
P (dv x ) = + g C y 

P (Cy) = “ dV X 
© 


p (dTy) = dVy 

P (dv ) =-gC x 
p (op = 1 dv y 
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Defining the state vector as 


dr 

dr 


x 


'y 


dv 

dv 


x 


the equations (2-8) and (2-9) can he written: 



"o 

0 

1 

0 

0 

0 ' 


"o 

o' 


0 

0 

0 

1 

0 

0 


0 

0 

• 

0 

0 

0 

0 

0 

g 


1 

0 

X = 

0 

0 

0 

0 

-g 

0 

X + 

0 

1 


0 

0 

0 

+ 1/R e 

0 

0 


0 

0 


_0 

0 

-1/R 

e 0 

0 

0 _ 


_0 

0 _ 


u 


u. 


(3-2-1) 


and for the measurement: 


N 


”l 

0 

0 

0 

0 

o' 

KJ 


L° 

1 

0 

0 

0 

oj 


x + 


x 

L v y. 


( 3 - 2 - 2 ) 


finally the matrices Q and R are : 



’n 

o' 


*R 

o' 

Q = 

_0 


R = 

_0 

R_ 


(3-2-3) 


Let us remark that here the first four rows and colum- 
ns of the matrix F are exactly the same as in model 1. This 
means that model 1 can he studied as a special case of mo- 
del 2 . This will he useful later on. 


3.4 Model 3 


The functional block diagram of the navigator is al- 
most the same as in the previous model, hut the misalign- 



2 



of H« I.M.5. 
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ment about the. z-axis is no longer neglected. Furthermore 
there is now an earth rate compensation which will intro- 
duce some other errors due to computed latitude feedback. 

„ 8 , 9 

See and . 

let us call "n" the navigation frame and '’cm" the con- 
trolled member frame - this is the platform axes -. With 
the same sign convention for the misalignment angles as be- 
fore (the angles are positive when the n-frame rotates about 

its positive axes to get aligned with the cm-frame) the di - 

cm 

rection cosines matrix can be written: 


1 

c 

-C 


1 

0 

o' 


0 

-C 



z 

y 







z 

y 

-c 

1 

c 


0 

r 

0 

— 

c 

0 

-c v 

z 


X 






z 


X 

c 

-C 

1 


0 

0 

1 


-C 

C v 

0 

L y 

X 







L y 

X 



or (T = I - C (3-3) 

where C is the skew- symmetric matrix associated with the 
rotation vector C (components • C x , C , C z ). Til is matrix is 
defined for any vector V by the relation : 


Yx W = [V] W 

Let us assume for the moment that position and velo- 
city are given in terms of latitude and longitude (angles 
instead of distances); thus: 

dv = p(L -L) fdr x = L q -L * (dv x )/p 

4 and 

dv y = P(l c -1) ldr y = l c -l = (dv y )/p 

In these relatione, the subscript "c" stands for "computed" 



22 


because 1 is one of the outputs of the I.R.S. and may dif- 

O 

fer from the true latitude L. 


A-l : the first assumption is to neglect Coriolis effects 
and vertical acceleration. 

Then the specific force in the n-frame is: 


f n = 


-R e p 2 1 
— R cosL p 2 l 

v 

g 


( 3 - 4 ) 


and the accelerometers output is the specific force coor- 
dinatized in the cm-frame which can he written 


f cm _ C cm f n 

Prom f n , vfi (angular velocity of the n-frame with res- 

' L vli 

pect to the earth frame, coordinatized in the nav. -frame) 
is given by the relation: 



cosl pi 


0 

-1/R e p 0 ' 


'-E e p 2 L 

(W n ) = 
v en'c 

-pi 

= 

1/R 

e P o 0 


-R cosl p 2 l 
e * 


-sinl pi 


0 

tanl /R p 0 

C G 


g 


or - (C’c = Ncf 1 (3 ~ 5 > 

In the expression of the performance function [pp] c , the com- 
puted latitude comes only into tanl . 

c 

But 

tan 1 = tan(L+dL) = tan 1 + dL/cos I +high order terms 

c 

We can, without introducing a great error, assume 
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tan 1 = tan I under the following assumption: 

u 

A-2 : the system is working far from the pole. 

At a latitude of 45 degrees, the resulting error is only 
2 dl (with dl in radians) and the approximation is. quite 
good. 

Prom , the velocity indication (in terms of Ion- 

gitude and latitude rates) is given hy: 


V 


pL 


0 -1 0 " 

_ v y_ 


Pi 


l/cos L c 0 0 


° r (I n ) 0 = [ p ] (3-6) 

The command angular velocity the space integrator 
must receive in order to operate properly is: 


iLa = o + 


en 7 c '-ie c 


(3-7) 


Thus (Wj e ) c = 


CJ . cos It 
le c 

0 


oo - „ sin L 
ie c 


must be added to W n 

— en 


The equation for the angular motion of the control- 
led member is now: 


—icm -in. •— ncm 


(3-8) 


where 


f cm 

■icm 


is the' commanded angular velocity given by 
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equation (3-7) and , which is the rate at which the 

controlled member goes away from the n-frame, is related 
to the correction angles by: 


— ncm 


w 

pc y 

.P c z . 


P C 


(3-9) 


Noting that wjjj = C° m wj n and using equations (3-5), 
(3-7) and (3-9), (3-7) can be written in the form: 


C 


cm 

n 


(W 11 + W 11 ) 

v ~xe — en' 


+ pC= [pf] 


c° m f n + 


(w n ) 

xe'c 


(3-10) 


Expanding sin L and cos L in terms of dL = 1 - L up 


to the first order', we get: 


(W? ) = W? - 

-ie c — ie 


co „• „ sin L 
xe 

0 

L “ie cosL . 


dL = W? + dW? 
— xe — xe 


Using (3-6), (3-10) becomes: 


[i - o] (wj e + v£ n ) *p£= [EF] [i - 0]f n . wj e + aw“ e 

(3-11) 

Some simplifications are possible. Indeed, from the prin- 
ciple of operation of the I.N.S. 

H £ n = Sen 


Further use of the relation [A]V = AxV = - YxA = - [v] A 
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yields finally: 


p c = awj e + [pf] f n sk 0 



(3-12) 


To be able to get a relation of the form (2-8), we 
need a first order differential equation between the seve- 
ral errors. Therefore we need to relate the terms 
[PP] [f n ] sk C to some other error terms. It turns out to 
be easily relatable to dv. Indeed, from (3-5) and (3-6) 


v + dv = [P] [PP] [I - C] f n 

dv = - [P] [PP] C f n 
dv = [p][PF][f n ] sk C 

We can develop this matrix product. Assumption A-2 allows 
us to write -l/cos L = -l/cos L (the error in- 

volved is of the same order as the one involved when wri- 
ting tan L = tan L) . 

C 

Furthermore, with the last assumption : 


Hence 

or 


A-3 : the vehicle is slowly moving on the surface of the 
earth 

the terms in pL and pi in the product [PP] ^f n j sk can be 
neglected. Thus 


f dv x" 


0 g/E e P 0 

1 


- g/R COS I p 0 0 

c 



(3-14) 
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Now we wish, to relate Jppj ^f^j slc C to dv. It is easy 
to see that 

0 cos L 
•1 0 
0 -sin L 

Substituting (3-14) into (3-12) yields finally: 





P (dr x ) 
P (dr y ) 



(3-16) 


Prom equations (3-13) , (3-15) and (3-16), defining 



where all the quantities involved 
are angles and angular rates 
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we can- write 


x = P x + G u 


with: 


P = 


0 

0 

0 

0 


■ cj. sin I 
ie 


0 1 
0 0 
0 0 
0 0 


0 

1 

0 

0 


0 

0 

0 


0 

0 


+ £_ 

e 


R cos L 
e 


0 0 cos L 


0 


0-1 


0 


0 


+ c*). sin L 

le 


0 


-u). sin L 
ie 


u. cos I 0 0 -sin L 
le 


and: 


0 


0 


-w cos L 
ie 


0 

0 

0 

0 

0 


+ (J. cos L 


0 


(3-17) 


G = 


0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

1 

0 

0 

0 


u = 


u 


x 


"y 


To get an equation of the same form as in the previous mo- 
dels, dr and dv must he expressed as functions of distance. 
Remarking that ■ 


'x 


L y J 


n. miles 


R 


0 


0 R cos 1 
e 

0 0 


0 


0 


0 

0 

R 

< 

0 


0 

0 

0 

R cos L 
e 


'x 


L7. 


minutes 


it follows that we can write : 


x = P x + G u 


with now 
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0 1 


0 0 


0 0 


0 0 


to . sin L 


0 0 + 


-"i e sln 1 


G> • cos L 


0 to. sin L 
le 

tan L 


(J. cos I 
le 


0 0 - 


-to. cos L 
le 


( 3 - 18 - 1 ) 

Then position is in nautical miles and velocity in nauti- 
cal miles per unit time. 

The measurement is m = H x + v with 


1 0 0 0 0 0 0 
0 1 0 0 0 0 0 


Finally, the matrices Q and R are the same as before: 


N 0 
0 N 


R 0 


0 R 


(3-18-3) 
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N and R are power spectral densities of white noises and 
thus are to he constant over some frequency interval. The- 
refore, the units must he: 

< 2 

(n. miles) .sec for R 

* 2 3 

(n. miles) / (sec) for N 

Let us recall the 3 assumptions made to derive this 
equation: 

A-l : neglect Coriolis effects and vertical acceleration 
A -2 : operation far from the pole 
A-3 : vehicle slowly moving. 
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CHAPTER 4 
FILTER EQUATIONS 


4.1 Introduction 

As it was pointed out in chapter 1, we are not inte- 
rested in finding the output of the filter for given I.N.S. 
and radar outputs, hut rather in measuring how accurate 
this indication is. Therefore, in this chapter, only the 
variance equation (the differential equation for P (t)) 
will he studied and solved. 

It follows from equations (3-1), (3-2) and (3-18) 
that the F matrices are very similar in all the cases. The 
12 equations governing the P (t) matrix in model 2 are a 

A 

particular case of the 28 equations of model 3; in the sa- 
me manner, the 6 equations of model 1 are a particular case 
of the 12 equations of model 2. 

Therefore, after an analytic solution for model 1, 
the only model that can he handled in a simple way, the 
equations for model 3 will he derived using 2 parameters 
"a” and "h” to allow a single study of the 3 cases. 
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4.2 Analytic solution of model 1 


The equations of this model are simple enough to he 
analytically solved. This allows us to find a closed form 
answer to the continuous filtering problem. 

let us go hack to equation (3-1-1). Since one of the 
main assumptions of this model is that there he no cross- 
coupling between the channels , the rank of all the matri- 
ces can he reduced by writing 


dr. 


dv. 


0 1 

0 0 

JF 


dr 


x 


dv. 


x 


0 

1 


u. 


G- u 


(4-1-1) 


and 


dr. 


x 


m - 


[x oj 


dr 

dr. 


x 


+ V. 


X 


H 


+ V 


(4-1-2) 


Thus the matrices Q and R become scalars N and R . 
If we let 


p x (t) = 


p ll (t) p X3 (t) 
P33OO 


the variance equa- 


tion can he written (dropping the variable t): 


1 

•■ci» • 
H 
H 

• 

p 13 


'0 

l” 


\l p l 3 ~ 


P 11 

p 13~ 


"0 0" 

■ 

p 13 

p 33_ 


0 

0 


_ p l 3 P 33 _ 

1 T 

- p 13 

p 33. 


1 °_ 


+ 

■* « 

0 

n["o 1 

-] - 

P 11 

p l3“ 


“l 

H 

O , 

» _ 
P 11 p 13 


1 

L 

J 

p 13 

p 33 


0 

R u J 

_ p 13 p 33_ 


(4-2) 
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4-2-1 Steady state 

The steady state solution is the solution of the 
right hand side of this equation when the left hand side 
is 0. ! 


The answer is straightforward 


'll ^ s . s . - 

VT R 4 

'4 

N 

'l3^s.s. 



'33^s.s. ” 

r‘ /4 



Thus, the steady state r.m.s. errors are: 


'in position RMX = 2^ R^N*^ 

^ in velocity RVS = 2*^ R ,/s N 3 / 8 (4-3) 


Now the optimum estimate is given hy equation (2-11) 


with 


K(t) = P(t) H* ( t ) R 


-1 


Let us call 

co„ - 

N 

‘4 

; then 

K 

3 a 

I 


n 

R 


s.s. 

<OnJ 


And: 


a *x 


0 1 


dr 

4. 


.<!>* 

i 


0 


dv 

_ - 


w n _ 


(m - [l o] 


dr 

dv 


(4-4) 


The optimum filter block diagram is given in figure 6 . 


4-2-2 Response to I.N.S. and radar noise 

Prom the diagram, it is easy to get the response of 


the filter to: 



ho X.N.S and 


rq daT noise & 


res p onse t. 
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#1 radar noise : 
on position 

on 'velocity 

#2 I.ff.S. noise 
on position 

on velocity 


e(£) 


p 1 

V u x 
n 

R 

p 1 ■ 

+ l/2“n 

P 

+ 

n 

e(v) 


"n p 



R 

_ 

P ■ 


P 

i 

+ “n 

e(r ) 


X 

p 



u/p 



P 

+ c»>* 

n 

e(v) 


1 +£w n 

P 

u/p 

= p*+ 


P 

n 


These response functions are plotted on figure 7 , 
4-2-3 Closed form solution for the free system 
The initial conditions can he taken as 


P x (-fc=0) = 


0 0 
0 0 

The variance equation for the free system is: 

• • 

P 11 p 13 

• • 

p 13 p 33 


or: I Pu = 2 P X3 

p 13 = p 33 
p 33 " H 

b 

With the initial condition 0, the solution to this equation 



"0 1 


p ll p l 3 ' 


~ p ll 

p 13 _ 

’o o" 


0 0 


0 0 


p 13 p 33 

+ 

p 13 

P 33. 

1 0 

+ 

0 H 
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is: 


P x (t) = N 


t 3 /3 


t a /2 


Then, "the r.m.s. errors are: 


on position RMX = 

on velocity RVX = 



t 


(H/3) ,y 2 t 3/jJ 
N l/l t & 


These results are plotted, on graphs 1 and 3 (upper curves) 
The steady state covariance matrix in continuous .fil- 
tering and the free system r.m.s. errors are about the on- 
ly things we can analytically study even in this simple 
model. As soon as we go into the transient solution for 
the continuous case, we get the following set of non-linear 
differential equations: 


P 11 ~ 2 p 13 Pll^ 
p 33 “ P 11 p 13^ 


p 13 

1*33 


N - P13/R 


These equations are not easy to solve (see 10 ) and 
better to solve them on a computer as a particular 
the most complicated model, as shown below. 


it is 
case of 


4.3 The variance equations 

Prom now on, the covariance matrix will be taken in 
the form: 


P 11 p 12 p 13 ' p 17 
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' (t)J with the usual definition 
for the state vector x(t) (chapter 3), the off diagonal 
terms are the correlations between the different state va- 
riables. 

For instance p-^g = E |^dv x C y j because all the sta- 
te variables are zero-meaned qualities (error terms). The 
diagonal terms are the variances (squares of the corres- 
ponding r.m.s. errors) for the same reason. 

For instance p^ = E [ dv yJ = (y-velocity rms error f 
let us come back to the system parameters matrix 
F of the 3 models, as they appear in equations (3-1-1), 
(3-2-1) and (3-18-1), and define three additional parame- 
ters a, b and c such that: 


Since P„(t) = E x(t) x 

A 


s 

a = 0 for model 1 and 1 otherwise 
4 b - 0 for model 2 and 1 otherwise 
c = ah 
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The differential equations for the covariance matrix will 
he derived with this expression for F . As for the number of 
equations necessary to study each of the cases, it appears 
to be : 


in model 1 - 6 equations corresponding to the 6 non zero 
quatities p^, p^, P 22 > ^24* P 44 p0r x- and 
y-channels. 


in model 2 - 6 more non zero terms: P^g? £36’ p 66 Por 
the x-channel and P 2 5> P 4 ^> P55 for the y-channel. There- 
fore, 12 equations are necessary. 


in model 3 - 28 equations because no quantity is a-priori 


zero. 



Equation (2-10) is: 

* 

P = F P + P F'+GQG' 

All the quantities have been defined and the result is the 
following set of 28 equations: 


P 11 = 2 p 13 

2 ) p 33 = 2ag p 36 + K 

■ 

3 ) P 13 = P33 + ag p lg 

4-) P22 = 2p 24 

5) P 44 = -2ag P 45 + N 

6 l £24 = £44 " ^ s _ p 25 

• .» S 1X1 Ii cl 

7) P 25 = P45 - t - P 3 2 + ~ bu ie 31n L p 26 

e E e 
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. w ie sin L 1 

8 ) P 45 = -g - b — g p 12 + — P 24 - b tJ ie sin 1 p 2g 

0 3 

. gj. sin L 2 

9) P 55 = -2b — g p 15 + ff- p 45 " 2bCJ le sinI, p 56 

© 0 

1 

1°) P 16 = P 35 - ff- P 13 + bcj ie sinl, p 15 + p 17 

o 

1 

11 ) p 36 = +6 p 66 - — P 33 + bu le sinl p 35+ bw ie oosL p 37 

0 

. 2 

12) P 66 = - E P 36 + 2bCJ le sinI p 56+ 2bw le oosL p 6? 

O 


1 3 ) P 12 = P23 + P14 • ' 

14 ) P 23 = P34 + s P 2 6 

15) p-]_4 = P 3 4 - g P]_5 

16 ) P 34 = S P 46 - g P35 

. to- sinL 1 

17) p 15 = P 35 ’P 1 i + r p l4“ U ie sin1, p l 6 

0 

. (j. sinL 1 

18) p 35 = gp 56 - — g P 13 + r P34- < * , ie slnI ' p 36 

0 0 

. 1 

19) P 26 = +P 46 - p 23 + o le sinI, p 25 +<^ ie cosl p 2? 

0 

1 

20) p 46 = -gp 56 - g P 34 + W 1 - e BinL p 45+ o. ie oosl p 47 

0 

. oj. sinL 1 1 

21) p 56 “ - -T ®16 + ff p 46- ^ie 31111 p 66" ff p 35 + CJ ie slnI ' p 55 

+ 4). e 00sL p 57 

. gj • cosL tanL 

22) p 1? » p 3? - -if- -P^- -ff— P 14 - ^ ie oosI p 16 

0 0 

. <o ie cosL tanL 

23) P 27 .= P 47 R" P 12 " ~R” P 24"" “ie cosL p 26 

0 0 

* ^ie 003 -^ lanL 

24 ) p 37 = gp 67 - — 5- — -p 13 - -ff— P 34 - ^ le =°si p 36 

e 0 

. co. cosL tanL 

"T 

25 ) p 4? = "gP 57 r- 'P14" "rT - ^- u ± e 00Bjl p 46 

0 0 
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. co- sinL 1 co. cost 

26) P 5? = £ P i7 + £ P 47 -Co ie smL p 6? p P 15 

e e 6 

tant 

p _ CO. COSt Pjt C 

R -^45 ie ^56 
0 

. 1 tO-: p cosL 

27) p 67 = - - p 37+ co ie sinL P 57 + co ie cosl p ?7 - — j -P 16 

e 

tant 

' p 46' t - J ie 00sI p 66* 

0 

co. cost tant 

28) p 77 = -2 ^ ?i 7 “ 2 “r P47” 2 ^ ie c °st pg 7 

0 0 

(4-6) 

In the second part of these equations (#7 through #12) the 
coefficient a has "been dropped; in the last part (#13 through 
#28) both coefficient a and b have been dropped. 

This system of linear first order differential equa- 
tions may be solved given some initial conditions and yields 
the r.m.s. errors in the state variable estimators. 


4.3.2 Continuous filtering compensation terms 


To get the equation (2-11) the quality 
M(t) = P x (t) H' R" 1 H P x (t) 


must be substracted from the foregoing equations. 

Of course M'(t) = M(t) so that only 28 terms must 
be computed as functions of the p. .. Assuming that M(t) 

J ' 


is written as 

M(t) 

tions is obtained: 



the following set of equa- 


1 ) m ll " (Pll + P 12 


2 ) m 


33 


= (p 13 + 


pJ 3 )/r 
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3) m-jj - (P 11 P 1 3+ Pl2 p 23^/ R ^ - 1 m 22~ ( p 12 + p 22 

5 ) » 44 = (p£ 4 + p| 4 )/r 6) ^ (p 12 p 14 + p 22 p 24 )/r 

7) ^2^ = (Pi2 p 15* ^ 22 ^ 2 5 m 45 = ^14^15"*" P 24 P 25) /// ' R 
9) m 55 = (p 45 + p| 5 )/R 10) m 16 = ^ p 11 p 1 6 + P 12 P 26 )/R 

U) nijg = ( p 13 p 16 + p 23 p 26 )/,R 12 ' m 66 = < - p 16 + P 26^ H 

13) m 2.2 = ^ p ll p 12 + P 12 P 22 )/^ ^4") m 23= (Pp2 p l3 + P 22 P ?3)^ R 

15) m 14 = (P 11 P 14 + Pi2 p 24 )/R 16) m 34 = (p 13 p 14 + p 23 p 24 )/R 

17) = ^ p nPi5 + p 12 p 2 5 ) ///r 1®) m 35 = ^ p 13 p 15 + P 23 P 25)^ R 

19) nigg = (P!2Pi6 + p 22 p 26) //R 20 ) m 46 = ^ P 14 P 16 + p 24 p 26) //R 

21) = (P 15 P 1 5+ P25 p 2 6^ 22 ) m 17 = ^ p ll p 17 + p 12 p 27^ 

23) ^2^ = ^ p i2 p l7 + p 22 p 27 ) //r ^ 4-) ^ P 13 P 17"*" P 23 P 27) ///R 

25) = ( Pq 4 Pq7 + p 24 p 27) ///r 2 ^) m 57~ ^ p 15 p 17"^ P 25 P 27) ///R 

27) m gy = ^ P 16 P 17 + p 26 p 27) //r 2 ^) m 77 = ^ P 17 + p 2 7 ) ^ R 

(4-7 ) 

These are the continuous filtering compensation terms writ- 
ten as functions of the p. ., some among them can he 0, de- 

cl 

pending on the model (the different models are delimited by 
the two dotted lines). 

To get the continuous filter variance equation, m. . 

ci 

is to be substracted from the right hand side of equation 
(4-6) : 

• • 

p ij “ ^ p ij)free system “ m ij (4-8) 

When the equation (4-8) are solved, the r.m.s. errors 


are given by: 
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'/» 11 = 

V? 22 = 

x/P^ = 


2? • Hi • S • 


If 


If 

n 


error in x-position = RMX 

" y-position = RMY 

v 

" x-velocity - RVX 

" y-velocity = RVY 


P 55> P 66 and 


P 77 


are 


the r.rn.s. 


mi s al i gnment angl e s . 


Let us find, in this case of continuous filtering, 
the form of the optimal filter. The optimum gains matrix 


K(t) = P 

(t) H' 

R" 1 

and is equal 

to : 


P 11 

P 12 


~ k ll 

k 12 


Pl2 

P 22 


k 12 

k 22 ' 

1 

p 13 

p 23 


k 13 

k 23 

K(t) = 

R 

p 14 

p 24 

= 

k 14 

k 24 


p 15 

p 2 5 


k 15 

k 25 


p 16 

p 26 


k 16 

k 26 


_ p 17 

p 27 


_ k 17 

ro 

-0 

i 


The optimum estimate obeys equation (2-11): 



9 

X - 

__ 

P x 

+ K 

( m - H x ) 

If “x 

is the measurement along 

the x-axis 

“V 

n 

U 

it 


,y-axis 


and if 

'V 

“x = 

m 

X 

- 

dr 

X 



■v = 

m 

y 

- 

...a 

dr y 


equation (2-11) can he written: 
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P(dr x ) 

P(dr y ) 

p (dT x ) 

p(dv v ) 


d ^x + k l A + k l2*y 

dv y + k 12 sr x + k 22 m y 

a gC y + k 13 m x + k 23 m y 
-ag C x + k 14 m x + k 24 m y 


p(C_) = -c 


siril &v ^ 

' R dr x~ a R ” c w ie sinL V k 15 m x +k 2 5 m ; 


P ( O y ) = 


p(3 z ) 


dv 


y 

/A 


y 


-a — - +c w. sinl C +c cosL C +k-, +k 0£r m 
R 0 le x le z 16 x 26 y 

tan! ^ 

~~ T " "■i 7 s 4 +k 27“y 


co. COSl ^ ^ ^ 

-C — dr x -c -g — dv y -c<o ±e cosL G y +k 17 E>k 9 , 7 m 


The signal flow diagram of the optimum filter for model 3 
is given in figure 12 , for purpose of simplicity any quan- 

'fci'ky k n(^ + k 2 _.m y (with k 24 = k^ 2 ) has been replaced 
by the number 11 j" . 

These quantities are computed from the 14 gains of 


equations (4-9) and the measured quantities m x and 
These computations are not shown on this figure. 


m 


4»3»3 Discrete filtering 

Between the measurements, during the operating time 
OPT, the covariance matrix obeys equations (4-6). At each 
measurement it is updated using equation (2-12); for this’ 
it is necessary to compute 


K = P 1 H* ( E P» H» + R ) -1 


It is easy to see that 
H P H» + R = 


p H+ E P 12 

p 12 e 22 + ^ 
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Since the general filtering theory assumes that R 
is a positive definite matrix and since the covariance ma- 
trix is non-negative definite, it is always possible to 
invert the matrix H P H' + R, even when the initial con- 
dition is P(t=0) = 0 , provided that we do not update P 
at the initial time. This is straightforward and can easi- 
ly be shown using the Schwartz inequality and the fact 
that P becomes more and more positive during the opera- 
ting time (from equation (2-10)). 

Denoting 


2 a 

D = det( H P H* + R) = P 11 P 2 2~ p i2 +R ( p ll +p 22) + R 


we can get : 


K(t) = - 


p ll^ p 22 + R ^ “ p 12 

Pl 2 (Pu+ R) 

- PuPx2 

P 12^ P 22 + R \“ p 12 p 22 

p 22^ p ll + 

~ p 12 

p l 3 ( p 22 + R) - P 12 P 23 

P 23 (PU + R) 

“ p 12 p 13 

p 14 ( ' P 22 + R ^ “ p 12 p 24 

P 24^ S 11 + R! 

~ p 12 p 14 

p 15^ p 22 + R ^ " p 12 p 25 

p 25 ( p ll+ H) 

“ p 12 p 15 

P 16^ P 22 + R ) “ p 12 p 26 

p 26 (p ll + E) 

” P 12 P 16 

p 17^ p 22 + “ p 12 p 27 

P 27^ P 11 + R ^ 

" p 12 p 17 


( 4 - 10 ) 


Once this is computed, the optimum estimate of x is: 


x = K 


“x 

“y 


The covariance matrix is updated by (I - KH) P-j 
in other words, is to he changed into : 



The new covariance matrix P(t*) "becomes the' new initial 
condition for the equations (4—6), until the next measure- 
ment is taken, when a new P(t^) becomes available. 

4.4 Computer program 

As it was already pointed out in part 4.2.3, the ana- 
lytic solution of these equations is not easy to obtain, 
even in the simplest model. Therefore, this work must be 
done on a computer. 

The computer program used for this paper is presented 
in Appendix B . 

The first program ■ (pages B-l through B-10) solves for 
the covariance matrix in the 3 modes : free system, conti- 
nuous, and discrete filtering; writes the r.m.s. errors, 
and punches some of these results for later use on a plot- 
ter. The main inputs are : -ALAT1 (latitude in degrees); 

TP (final time); AQ (inertial navigator noise in ft /sec ) 
AR (radar noise in ft .sec), and three different operating 
times 0PT1, 0PT2, OPT3 in sec. 
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CHAPTER 5 


RESULTS 


The computer program used to get the results of this 
chapter is given in Appendix B . It was worked out on an 
I.B.M.-360 in the M.I.T. Computation Center.' 

Por the discrete filtering, 3 operating times were 
studied: 30 seconds, 3 minutes and 18 minutes. 

Por any of the 3 modes, several cases for both the 

I.'N.S. and the radar noise have been worked out; N, the I. 

.2 l . .2 . .3 

H.S. noise, ranging from 10 to 10 (feet) /(sec) (po- 
wer spectral density of white noise over some frequency 

+ 6 

range) and R, the radar noise, ranging from 10 ' to 10 

2 

(feet) .sec. These values seem to cover all the practical 
cases . 

Before exposing the results, it is better to give 
some explanations about the curves of Appendix A . * 

The first 10 graphs are actual outputs from the com- 
puter and represent the variation of RHX (r.m.s. error in 
x-position indication) and RVX (r.m.s. error in x-velocity 
indication) as functions of time in minutes. The units are 
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nautical miles for RMX and feet/sec for RVX. . 

Graphs 1 and 2 (3 and 4) are RMX (RVX) for the free 

i 

system and for the 3 models. 

Graphs 5 through 12 give RMX and RVX for models 1 
and 2 only. There are 4 curves on each gx*aph: 1 for the 
continuous filtering and 3 for each operating, time in the 
discrete filtering. 

Graphs 13 through 19 have been established from the 
printed outputs 'and the units are : feet for RMX and ft/see 
for RVX. Graphs 13 through 17 are concerned with continuous 
filtering only and chart 17 yields the expected RMX and 
RVX for given I.R.S. noise N and radar noise R, both in u- 
sual units . 

The last 2 graphs are concerned with discrete filte- 
ring. from them we can find the influence of the operating 
time on both RMX and RVX given some values for the r.m.s. 
errors in the continuous filter (see below). 

5.1 free system ( graphs 1, 2, 3 and 4) 

According to the results of part 4.2.3> RMX and RVX 

3/, I/, 

increase as t for RMX and t for RVX for a given 
inertial noise in model 1 . This is observable on graphs 1 
and 2 -for RMX, 3 and 4 for RVX, where the upper curve re- 
presents model 1. 

The 2 other curves on these graphs represent the free 
system r.m.s. errors and are almost the same. The improve- 
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ment in the error by comparison with the model 1 is due to 
the Sculer oscillation (indeed the time "between 2 consecu- 
tive points of both curves where the slope is minimum ap- 
pears to be 42 minutes, half of the Sculer period). This 
improvement with respect to model 1 reaches 80 percent for 
RMX an d 30 percent for RVX after 84 minutes. 

The tiny difference between the curves representing 
models 2 and 3 was also observed for other values of N and 
clearly shows that model 3 is not better than model 2 when 
there is no filtering. 

Furthermore, up to 14 minutes for RMX and 8 minutes 
for RVX, the 3 curves are almost coincident. This means 
that the 3 models yield the same r.m.s. errors in free mode 
up to 14 minutes for position error and 8 minutes for ve- 
locity error . 

Of course, since no radar is used here, the graphs 
do not depend on R. Furthermore, looking at the upper cur- 
ve on graphs 1 and 2 as well as 3 and 4, the influence of 
N can be found to match equation of part 4*2.3 • In graphs 
1 and 2 for instance, the scale of the RMX-axis is multi- 
plied by 10 while the noise is multiplied by 100. The same 
thing stands for RVX. And since the curves have exactly the 
same shape in both cases, the r.m.s. position and velocity 
errors are proportional to at a given time . 

The general result for the free system can be stated 


in the following way; 
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for model 1, the r.m.s. errors are given by 
RMX = (N/3) Vi t 3/i 

Vi *!x (5-D 

RVX = N Z t Z 

k 

At any time, there is no difference between model 2 and 3. 
The 3 models are equivalent up to 14 minutes for position 
error and 8 minutes for velocity error. 

5.2 Continuous filtering 

The curve representing this mode is the lowest one 
on graphs 5 through 12 . It is the only straight line- on 
all these figures. Because the steady state error in conti- 
nuous mode is very small and not easily readable on. the 
plots, the results are given on page A-4, with the preci- 
sion obtained on the computer. Bor all the cases that have 
been studied, the 3 .values of RMX (in nautical mile) cor- 
responding to the 3 models and the 3 values of RVX (in 
feet/sec) corresponding to the 3 models are shown. 

The largest difference between the 3 models appears 
to be .04 percent and is indistinguishable from the trun- 
cation errors in the results. 

Therefore, and this conclusion is an important one, 
the steady state position and velocity r.m.s. errors do not 
depend on the model chosen to represent the I.N.S. in the 
case of continuous filterin. 
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It could "be expected from part 5.1 (free system) , 
where models 2 and 3 appeared to yield the same error, that 
there were no difference between the 2 models. The important 
fact is now that model 1 yields also the same result. 

The steady state errors for model 1 were shown in 
part 4.2.1 to obey equations (4-3) • To argue the model inde- 
pendence, the same formula can be found from graphs 13 - 16 
where the results for model 2 are plotted. 

Graphs 13 and 14 show the -variation of RMX with R for 
different values of R and with R for different values of R. 
On logarithmic paper, the curves appear to be straight pa- 
rallel lines. Measuring the slopes yields the following 
dependence formula; 


1 3 

log RMX = — log R + — log R + d 

8 8 

1 

and the constant d turns out to be — log 2 

4 


Thus : 


RMX = 2 V ‘/V /S 


A similar analysis for RYX. yields: 
IrvX = 2^ R^ 8 R^ 8 


In these equations the units are: 


RMX in feet ; RTX in feet/sec 

2/3 a 

R in feet /sec ; R in feet .sec 


(5-2) 

(5-3) 


These results could also have been found by a dimensional 
analysis, R and R beeing the only parameters that can in- 
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fluence the errors. 

The last thing to be said about this continuous fil- 

i 

tering is that, as it was expected from figure 6 in model 
1, the time necessary to reach the steady state is small 
and does not exceed 5 minutes in the worst case presented 
on graph 3 • 

Graph 15 gives the position and velocity r.m.s. er- 
rors in. feet and feet/sec for given values of N (I.N.S. 
noise) and R (radar noise). 

5»3 Discrete filtering 

Let us recall that this case was studied with 3 dif- 
ferent operating times: OPT = 30 seconds, 3 minutes and 18 
minutes, which seem to cover a good part of the permissible 
range . 

The corresponding curves are the 3 upper curves on 
graphs 5 through 12, and it is obvious that the larger the 
operating time is, the higher the corresponding curve goes. 

Let us note that, especially for the lower two OPT, 
the curves are not quite representative of what happens. 
Indeed we have assumed that updating after each measurement 
was an instantaneous operation, so that the curve should 
drop with an infinite slope at each measurement time. 

This is not the case because, in order to save some com- 
puter time, we limited ourselves to 3 output values per 
operating time, and the 3 values could not be chosen to be 
just prior to and just after the measurement. 



52 


Furthermore, it is very difficult to analyse these 
results because the "mean value" for each case is uneasy 
to obtain, either from the graphs or from the printed out- 
puts . 

Therefore, let us see if it is possible to relate 
any discrete filtering problem to what will be called its 
" continuous approximation " . 

Let us start with a discrete measurement process 
obeying the usual equations i 

x(t) = F(t) x(t) + Gr(t) u(t) 

[m(t n ) = H(t n ) X(t n ) + v(t n ) 

where u(t) is a white noise and v(t ) an independent Gaus - 
sian process such that: 

'e jju( t ) u’(t+s)]= Q(t) S (s) 

4 E[v(t n ) v‘(t n )J = V = constant 
E j^v(tn) v‘(t n+ s)] = 0 if s is not 0. 

Since the noises are zero-meaned, V is. the variance of the 
radar noise. 

This discrete measurement process can be approxima- 
ted with an equivalent continuous one defined by 

m(t) = H(t) x(t) + v(t) 

where now v(t) is a white noise obeying: 

EjV(t) v'(t+s)] = R $ (s) with R constant. 
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Furthermore let us assume 

R = V . OPT (5-4) 

The units are right since V is a variance and R a power spec- 
tral density for a white noise. This really means that if 
the measurements are taken twice faster, the measurement 
error covariance must "be twice larger to be approximated 
by the same continuous process. 

Then it can be shown that, provided that the operating 
time is not too large , any discrete measurement process 
can be approached by a continuous one which appears as the 
mean of the previous one. This can be understood by deri- 
ving the continuous case variance equation from the discre- 
te case scheme. 

Between two measurements times t n _-^ and t , equation 
(2-10) can be written with the state transition matrix 
<$> (t^tj) as: 

p(t~) = <3>(t ,t ) p(t -, ) <£>' (t ,t ) 
v n' n* n-1' ' n-1' ^ v n’ n-1' 

+ I t>(t n ,s) S(s) Q(s) G'(s) <J>‘(t n ,a) as 

^n-l (5-5) 

If and P are continuous, Taylor expansion yields: 

<J) (t n ,t i) = I + P(t n _ 1 ) At + higher order terms in At 
At is here the operating time \ L ”' t n-l* 


Then, using the mean value theorem with 
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(5-5) becomes: 


= p <Vf + 5 <Vi) p <Vi> ^ + «Vi> 


+ G(z) Q(z) G'(z) At + higher order terms 


(5-6) 


At measurement time we update P(t ) through 


p(tp = p(t n ) - p(t;j h' 


n 


n' 


R 

H P(t ) H' + - 
n At 


-1 


■H P(t") 


This yields finally the equation: 

P(t n )' = p(Vi> + [ p< ' iil1 - 

+ G(z) Q(z) G'(z) - P(t") H*[h P(t“)H'At + vj H P(t~^ 

(5-7) 

And in the limit when At approaches 0, we get 


n ) P(t 1 ) + P(t , ) P» (t -, ) 
1' v n-1 v n-1 7 n-l / 

-1 


At 


P(t) = F P + P F’ + G Q G’ —PH 1 R 1 H P 


which is the continuous filtering variance equation. 

Thus, under assumption (5-4), any discrete problem 
can be approximated by a continuous one since has the 
same effect as 7. It happens that, for not too large an 
OPT, this continuous approximation looks as the mean of 
the discrete process and is therefore the best measure of 
its accuracy . 

Prom our point of view, the accuracy of the discrete 
filtering system with operating time OPT and error variance 
V is the same as the accuracy of a continuous filtering 
system with noise power spectral density R = V . OPT 
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Let us see now the significance of the limitation 
on the time "between the measurements. Two reasons can "be 

i 

thought of : 

1) if the discrete process yields the same results as the 
continuous one, this means that the r.m.s. errors do not 
depend on the model in the discrete process (since this 
was shown in the continuous case). But it was pointed out 
in part 5.1 that the 3 models are completely equivalent 
in free mode only during a time lower than 14 minutes for 
position error and 8 minutes- for velocity error. Since this 
free mode is precisely used "between the resets, it seems 
consistent to take as time limits for this continuous app- 
roximation of the discrete process the 2 numbers: 14 minu- 
tes for position and 8 minutes for velocity. 

2) what is really implied by the relation (5-4) is to re- 
place,' on a plot of the autocorrelation function of the 
radar noise as a function of time, an impulse (rectangular 
shape) of area R by a triangular curve of hight 7 and area 
7 . OPT = R . Although these two curves are not the same, 
they can yield the same final result if the system is una- 
ble to distinguish between both shapes. And this happens . 
if the spread of the autocorrelation function is small 

by comparison with the time constant of the system. The 
practical limitations appear to be the same as before: 14 
minutes for position and 8 minutes for velocity, as long 
as the Sculer oscillations do not change the response too 
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much. 

Under these limits, from (5-2), the steady state er- 
ror of the continuous process does not depend on the model. 
Thus the accuracy of a discrete process does not depend on 
the model . This, of course, stands as long as we can speak 
of a mean for this process; that is to say that the opera- 
ting time must be small by comparison with the mission ti- 
me . 

Graphs 18 and 19 give the r.m.s. errors in position 
and velocity for discrete processes of given operating ti- 
me and steady state errors of the corresponding ( not equi - 
valent ) continuous process. They are to be used- in con- 
junction with graph 17 in the following way: 

for any I.N.S. and radar noise powers N and R, graph 
17 gives the resulting position and velocity r.m.s. errors 
in the steady state. With these values of RMX and RVX, 
graphs 18 and 19 give the corresponding r.m.s. errors of 
a discrete filter of given operating time OPT (in seconds') 
using the same I.N.S. and radar. 

The accuracy of the continuous approximation can be 
checked on graphs 6 and 11 where the mean values of the 30 
seconds, 3 minutes and 18 minutes discrete processes have 
been plotted. It is to be remembered that the shapes of the 
curves are not exact and the error increases with the ope- 
rating time. 

Anyway, it appears that the approximations are quite 
good for the 30 seconds and 3 minutes cases. As for the 18 
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minutes process, while the position error curve is still 

acceptable, this is no longer true for the velocity error 

\ 

curve. This fact is in accordance with the operating time 
limitations previously introduced. It can he found conve- 
nient to consider 10 minutes as the limit on the opera- 
ting time for hoth position and velocity informations. 

5.4 Summary of the results ; Conclusion 

1. For a use of the inertial navigator alone, without ex- 
ternal position information, the 3 models are equivalent 
during a time not exceeding 10 minutes when they start 
from the same perfect initial state. After 10 minutes, the 
Sculer oscillations attenuate the errors in models 2 and 3 
which are always equivalent. The errors in model 1 are gi- 
ven hy equations (5-1). 

2. When use of a position information is possible in a 
continuous way, the model chosen in the Kalman filter to 
represent the I.N.S. is of no importance. This means that 
a simple model consisting of 2 accelerometers kept rough- 
ly aligned with the north and east axes yields the same 
accuracy as a more sophisticated one (hut becomes very had 
if the external information happens to he lost) . 

For any model, the steady state r.m.s. errors are 
given hy equations (5-2) and (5-3). 
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3. As long as the correlation time of the external device 
(radar) is small by comparison with the time constant of 
the system, any discrete measurement process can be appro- 
ximated by a continuous one in the sense that the continuous 
process appears as the mean of the discrete one. 

As long as the operating time is smaller than the 
time necessary for the Sculer oscillations to attenuate 
the errors, the 3 models yield the same final error. 

Since this last constraint is less drastic than the 
first one, the result can be, stated in the following way: 
as long as the operating time does not exceed 10 minutes, 
the position and velo'city' r.m.s. errors in a discrete 
measurement process do not depend on the model and can be 
approximated by the errors in a continuous scheme related 
to. the discrete one by the equation (5-4). 

2,3 

The final equations are, with N in feet /sec , R 
in feet* , OPT in seconds , RMX in feet and RVX in 
feet/sec : 
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APPENDIX A 


GRAPHS 


The units in the following graphs are the following: 

(f eet)/seo'. for velocity error RVX 

nautical miles in graphs 1 through 12 as well as on 
page A-3 and feet otherwise for position error RlflX. 

On graphs 1 through 12 which are output from the com- 
puter, some indications appear inside the hox in upper left 
corner: the mode (1 for free system and 3 for discrete fil- 
tering); whenever the model does not appear, 3 curves cor- 
responding to the 3 models are plotted. The other values 
are N in (feet) / (sec) , R in (feet) .sec and 
the 3 operating times in minutes. 
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APPENDIX B : COMPUTER PROGRAM 


PROGRAM FOR COMPUTING POSITION AND VELOCITY R.M.S. 
ACDITTONAL OUTPUT : MISALIGNMENT ANGLES R.M.S. 


DIMENSION TI(50) ,X)l( 5C ) , X2 1 1 50 » , X 3 1 ( 50 ) ,X12(50) , X 22( 50 ) , X32 < 50 ) ,V 
lll( 50) ,V2l( SO) , V31I50) ,V12( 50) ,V2?( 50),V32(50) ,X131 (50) ,X132(50 ) ,X 
2133(50 »X231(5C)*X2 32(50) »X233( 50 ),V131( 50 ),V132(50)»V133(50),V231 
3(50), V232(50) ,V233(50) 

COMMON A, R,Y ( 28) , YP( 28 » , R, N, J,NS 

COMMON G,RE,0X,r7,0RX,ORZ,TR,Q,M0nF,PAS,E' , SI,SEUIL( 23) ,T 
COMMON CK( 14) 

******************************* 

READ AND WR ITF DATA 

READ (5,1^01) ALAT1 ,EPSI , P A 1 , TF , DT W , SE , NPU , TPU1 , K , MODEM I , MODEM A 

1001 FORMAT (F5. 1,E8.3,F8.3,E8.3,F8.3,E8.3,12,E8.3,I3,I2,!2) 

RE=6366000. 

G=9.8067 

01E = 3. 1416/43200. 

FC1=3. 28084 
FC 2=1 852. 

Al AT=( 3. 141 6/180. )*AL ATI 

1 READ (5,1002) AO , A R , OPT 1 , 0PT2 , OP T3 

1002 FORMAT ( 5E8. 3 ) 

K=K* 1 

WRITE (6,1003) ALAT1,EPSI,PA1,TF,0TW 

1003 FORMAT ( 1HI , 1 2X , 14H0PT IMUM Ml X I NG// / 1 2X , 1 1HL AT ITUDF = F 5 . 1 , 1 2 X, 7H EP 
1SI = E 1 1 • 4 » 1 2 X , 6HP A S = E11.4,12X, 5HTF = FI 1 .4 // 20X , 6HDT W = Ell. 41 

IF (ALAT.E0.90. ) GO TO 22 
***************************** 

COMPUTE USEFUL TERMS 

ALAT= (3. 1416/180. ) *AL AT 
OX=OI E*COS( ALAT ) 

0Z = -0!E*SIN( ALAT) 

ORX=OX/RE 

ORZ=OZ/RE 

TR* SI N ( ALAT )/ (R E*COS( ALAT ) ) 

Q= AQ/ ( FC1 *FC1 ) 

R=AR/ ( FC1*FC1 ) 

******************************* 

BEGIN EACH CASF ; WRITE INITIAL CONDITIONS 

DO 35 MODE L = 1 » 3 

DO 21 MOD E= MO OEM I , MODEM A 

0PT=0PT1 

2 P AS =t> A 1 
A= 1 • 

8=1. 

TW=DTW 

TM=OPT 

T PU=T PU1 
T =0 . 

WRITE (6,1004) AQ, AR, OPT , MODEL, MODE 

1004 FORMAT ( 1H1//12X, 4HN = E 1 1 . 4 , 14H ( FEET ) 2/ ( SEC » 3 , 1 2 X , 4HR = Ell.4,llH( 
1FEET) 2*SEC , 1 2X, 6H0PT = E11.4,4H SEC //20X , 8HM0DEL = I 2 ,20X, 7HM0DE = 
2 12) 
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B. 


GO TO (3,4,5) , MODEL 

3 N=6 
A=0 . 

GO TO 6 

4 N=12 

R=0. 

GO TO 6 
* 5 N=28 

. 6 DO 7 1=1,28 
Y ( I ) =0 . 

7 SEUIL < I ) = S E 

GO TO (9,10,11) , MODEL 
" 9 WRITE (6,1005) 

LOO 5 FORMAT(/20X,4HT = ,20X,6HRMX * ,2^X,6HRVX = //) 

WRITE (6,1009) T,Y(1),Y(2) 

GO TO 12 

10 WRITE (6,1006) 

L006 FORMAT ( 10X, 4HT = ,20X,6HRMX = ,20X,6HRVX = ,20X,6HRCY = /) 

WRITE (6,1010) T,Y( 1 ) ,Y(5) ,Y(12) 

GO TO 12 

11 WRITE (6,1007) 

L 007 FORMA T(6X, 4HT = ,10X,6HRMX » ,8X,6HRMY = ,8X,6HRVX * ,8X,6HRVY * , 
18X,6HRCX = , 8X , 6HRC Y = ,8X,6HRCZ = /) 

WRITE (6,1008) T,Y(1),Y{4) , Y( 2 ) , Y( 5 ) , Y ( 9 ) , Y( 12) ,Y(28) 

jfc jjcjjc j(c j$cjJc#j$c 5 $: ;$t3!c;6c3t3}c:jc:(c 

CALL DIFF. EQUA. FOR INITIAL DERIV, 


12 M=0 

CALL DAUX 

** ****** * * ** * *** * ********* * * * * * 

CALL SUBROUTINE OF INTEGRATION 

13 CALL KUTAM- 

IF (M.GE.50) GO TO 210 

* 3}! $ * * £ * * * * * * sfr * * * * * * * * * * * * * * * * * * 

IF TEST SI CALLEO TOO MANY TIMES »NEXT CASE 

IF (NS.GE.50) GO TO 1 

* * * * ** * ** * * * * * * * * * * * * * * * * * * * * ** 

IF INTEGRATION STEP SIZE TOO L ARGE , INI TI AL VALUE 

IF (PAS. GT. PAD PAS=PAl 

aft * * * * * * * * * 5ic * * * * * * * * * * * * * * * * * * * * 

TEST ON FINAL TIME 
IF (T-TF) 14,210,210 

* * * * * 3fc * * 4c * * 3$c i}e * * * * * * * $ * * * * * * * * * * 

TEST ON' TIME AND MODE TO CALL DISCRETE CASE COMPENSATION 

14 IF ( T-TM.GE . 0.*AND.MODE.EQ.3) GO TO 36 
GO TO 16 

36 CALL UPDAT 

IF (J.GE.5) GO TO 23 
TM=TM+OPT 

********************** ********* 

AUX. CALCULUS {WRITE RESULTS 
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16 RMX=SQRT(Y(1))/FC2 
RVX=SQRT(Y(2) >*FC1 
RMY=SQRT(Y(4))/FC2 
RVY=SQRT (Y(5))*FCl 
RCX=SQRT( Y I Q ) )*1000. 

RCY=SQRT(Y( 12) ) *1000. 

RCZ=SQRT( ABS(Y(28> I >*1000. 

a****************************** 

IF NPU=1 AND T.GE.TPU , PREPARE PUNCH 

IF (NPU.NE. II GO TO 33 
IF (T.LT.TPU) GO TO 33 
M=M+1 

TPU=TPU+TPU1 

TI(M)=T 

GO TO (40,50,60) , MODEL 

40 GO TO (41, 42,43) , MODE 

41 XI 1 ( M ) = RMX 

V 1 1 ( M ) =RVX 
GO TO 33 ■ 

. 42 X 1 2 ( MJ^RMX . 

V12 (M ) =RVX 
GO TO 33 

43 IF (DPT— 0PT2) 431,432,433 

431 X131(M)=RMX 
V13K M)=RVX 
GO TO 33 

432 X 132 ( M )=R MX 
V132( M)=RVX 
GO TO 33 

433 XI 33( M)=RMX 
V133 ( M) =RVX 
GO TO 33 

50 GO TO (51, 52, 53), MODE 

51 X2 1 ( M ) = RMX 
V21 ( M ) — RV X 
GO TO 33 

52 X?2(M)=RMX 
V22 (M ) =RVX 
GO TO 33 

53 IF (OPT-OPT?) 531,532,533 

531 X231 ( M) =RMX 
V231( M ) =RV X 
GO TO 33 

532 X?32( M ) = RMX 
V.? 32( M)=RVX 
GO TO 33 

533 X233(M)=RMX 
V233 ( M )=RVX 
GO TO 33 

60 GO TO (61,62,33) .MODE 

61 X31(M)=RMX 
V31 (M)=RVX 
GO TO 33 

62 X32 (M) =RMX 
V32(M)=RVX 

******************************* 

IF T.GE.TW , PRINT RESULTS 
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33 IF(T.LT.TW) GO TO 13 
GO TO ( 18 T 19,34 ) , MODEL 

34 WRITE ( 6 » 1008 ) T , RMX , RMY , RVX , RVY , RCX , RCY T RC l 

1008 F0RMAT(4X,E11.4,3X,E11.4,3X,E11.4, 3X , E 1 l . 4, 3X , E 1 1 . 4, 3X , E 1 1 . 4, 3X , 
1EI1.4,3X,EU.4> 

TW=TW+DTW 

GO TO 13 . 

ur 18 IF (RMX. NE. RMY. OR. RVX. NE . RVY) GO TO 20 
A WRITE (6,1009) T, RMX, RVX 

1009 F0RM4T(18X,E11.4,13X,E11.4,15X,E11.4) 

TW=TW+DTW 

GO TO 13 

U “ 19 IF (RMX. NE. RMY. OR. RVX. NE. RVY) GO TO 20 
WRITE 16,1010) T,RMX, RVX,RCY 

1010 FORMAT! 8X f E 1 1 .4, 13X ,E1 1 . 4 , 15X, El 1 . 4, 15X,E11.4) 

TW=TW+DTW 

GO TO 13 

20 WRITE (6,1011) 

1011 FORMAT (SOX, 24HCHANNELS NOT INDEPENDENT) 

210 WRITE (6,2100) M 
2100 FORMAT ( 2X,4HM = 13) 

IF (MODE-2) 21,21,214 
214 IF (0PT-0PT2) 211, 212, 213 

211 OPT=OPT 2 
GO TO 2 

212 OPT^OPT 3 
GO TO 2 

213 GO TO 35 

21 CONTINUE 

35 CONTINUE 

# if if if & $ sjt * * # >$: if # # if # if # # 

IF NPU = 1 , NORMALIZE RESULTS OF EACH STEP 

IF (NPU.NE.l) GO TO 71 
XM1=0 • 

XM2=0 . 

VM1 =0 . 

VM2=0. 

XM13=0. 

XM23=0. 

VM13=0. 

VM23=0 . 

DO 100 M=1 , 50 

EXl^AMAXKXll(M) ,X21(M) ,X31(M) ) 

EX2=4MAX1(X12(M) ,X22( M) , X32 ( M) ) 

EV1=AMAX1(V11(M) ,V21(M) ,V31(M) ) 

EV2=AMAX1(V12(M),V22(M),V32(M) ) 

EX13=AMAX1(X131(M),X132(H) ,X133(M) ) 

EX23=AMAX1( X231 ( M ) , X232 l M ) , X23 3 ( Ml ) 

EV13=AMAXI (V131(M) , V 132 ( M ) , V133 ( M ) ) 

EV23=AMAX1(V231(M) ,V232 (Ml , V233(M) ) 

IF (EX1.GT.XM1) XMl-EXl 
IF (EX2.GT.XM2) XM2=EX2 
IF (EV1.GT.VM1) VM1=EV1 
IF CEV2.GT.VM2I VM2=EV2 
IF (EX13.GT.XM13) XM13=EX13 
IF (EX23.GT.XM23) XM23=EX23 



78 


B-5 


IF (EV13.GT.VM13) VM13=EVI3 
IF (EV2^.GT.VM23) VM23=EV?3 

100 CONTINUE 

D n ioi M= 1 , SO 
XI 1(M)=X11 
X?l(M)=X?l{M)/X M l 
X 3 1 ( M)=X31 ( M ) /XM1 
X12(M)=X1?(M)/XM2 

X??(M)=X2? IN) /XM? 

X32(M)=X32(M) /XM2 
Vll (M)=VI 1 ( M ) / VM 1 
V21(M)=V2H M )/VM1 

V3 1 ( M ) = V3 1 ( M ) /VM l 
V12(M ) =V12 I M I /V M2 
V2?(M)=V2?(M)/VM? 

V32(M)=V32(M)/VM2 
XI 311 M I =X I 31 ( M) / XM 13 
X132I M ) =X 132 ( M) /XM 1 3 
X 1 331 M )=X 133 1 M ) /XM1 3 
X231 ( M)=X?31 CM) /XM23 
X232I M)=X232(M) /XM23 
X233 { M ) =X233I M) /XM23 
V13H M ) = V 1 3 1 IM) /VM13 
V132I M) = V1 32 ( M ) /VM 1 3 
V133I M )=V133 CM) /VM13 
V?3l(M)=V231{M)/VM23 
V232( M)=V232(M) /VM23 
V?33{ M) =V?33{ Ml /VM23 

101 continue 

*3 fle*# £ >}: * s}: 3|{ 3|t 

IF NPU = 1 , PUNCH RESULTS OF EACH STEP AND NORMALIZATION CONSTANTS 
PUNCH 1C20,K,AQ,AR 

1020 FORMAT 1 5HCASE I?,2X,4HN = Ell.4,3X,4HR = Ell. 4) 

PUNCH 1021, 0PT1 ,0PT2,0PT3 

1021 FORMAT (7H0PT1 = El 1 . 4 , 3 X ,7 HOP T2 = E 1 1 . 4, 3X ,7H0P 13 = EH. 4) 

PUNCH 1022,XM1, XM2,VMl, VM2 

1022 <=ORMAT I6HXM1 = E11.4,6HXM2 = E11.4,6HVM1 = E11.4,6HVM2 = Ell. 4} 

PUNCH 1023 ,XM13,XM?3» VM13,VM23 

1023 FORMAT (7HXM13 •= E11.4,7HXM23 = E11.4,7HVM13 = E11.4,7HVM23 = Ell. 4 
1 ) 

nn 102 m= i ,50 

PUNCH 1024, TI { M) ,X11 (M),X21<M) ,X31{M) , X12( M ) , X2 2 { M ) , X32 ( M ) ,V11( M) , 

IV? 1 (Ml ,V31( M) ,V12(M) , V22 (M) ,V32< M) ,K,M 
PUNCH 1024,111 M) ,X1 31 I MI , X132 ( M I , X133{ M ) , X23HM) , X232( M) ,X233( M> , V 
1131 (MI ,V132(M) ,V133(M) ,V23l ( M) , V232( M) ,V233( MI , K, M 

1024 FORMAT (F6.C, 12FS.3, 6X, I2,2X,I3) 

102 CONTINUE 
71 GO TO l 

2? WRITE (6,1012) 

101? FORMAT (25X,20HLATITUDE NOT ALLOWED) 

23 WRITE (6,1013) MODEL, MODE 

1013 FORMAT ( 10X, 7HM0DEL ,I?,7H MODE , I 2 , 10HI M POSS I BL E ) 

STOP 

END 
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SUBRCL'TINE CAUX 
DIMENSICN C T ( 2 8 ) 

CCPPCN A,E»Y(28),YP(26)»R»N*JtNS 

CCNNCN G,RE,CX,CZ t CRX , C R 7 ,T R , C » MCC E , PAS T E PS I , S EU IL { 28 ) , T 
COPMZN CK { 1 4 ) 

C=A*E 

FREE SYSTEM EQLA. 

YPU) = 2.*Y (3) 

Y P ( 2 ) = 2 • * A*G * Y ( llt+O 
YP(3 )=Y (2)+A*G*Y(10) 

YP ( 4 ) =2 . * Y < 6 ) 

YP(5) = -2.*A*C-*Y(8)+Q 
YP ( 6 ) =Y (5 )-/ ^G*Y (7 ) 

IF (A.EO.C.I GC TC 1CI 

Y P ( 7 ) = Y (8 )+fi*DRZ*Y { 13 ) + Y ( 6 )/RE+B*GZ*Y( IS) 

YPm=-G*Y(9) + B*CPZ*Y (15 ) + Y (5)/RE+B#CZ*Y{20) 

YP(S) = 2.*e*0RZ*Y{ 17)+2.*Y< 8 ) /RE + 2 .*B*CZ*Y ( 2 1 ) 

YP(10 )=Y<11 1-Y(3 )/RE-B*OZ*Y( 17 }+B *OX*Y ( 2 2 ) 

YP( 11 )=G*Y( 12 )— Y{ 2)/RE-E*CZ*Y{ 18 ) + E*CX*Y (24 ) 

Y P ( 1 2 ) =— 2 .*Y( 11)/RE-2.*B*QZ*Y( 2 1 ) +2 ,*B*G X* Y ( 2 7 ) 

IF (E.EC.O.) GC TC 101 
YP ( 1 3 ) = Y ( 1 4 ) + Y ( 15) 

Y P { 14 )=Y ( 16 )+G*Y{ IS) 

YF(15)=Y( 16 ) -G# Y (17) 

YP ( 16) = G*Y(2C) -G * Y ( 1 8 ) 

YP{ 17 )=Y { 18 )+ORZ*Y { l ) +Y { 15 )/RE+OZ*Y{ 10) 

YP< 1E)=G*>( 21>+CRZ*Y (3 ) + Y ( 16 ) /RE+ C Z*Y ( 1 1 ) 

YP ( 19 ) = Y ( 20 )— Y ( 14)/RE-0Z*Y( 7>+GX*Y(23) 

YP(2C)=-G*Y (21 )-Y (16 )/RE~CZ*Y (8 )+GX*Y{2 5.) 

YP( 2 1 )=ORZ*Y ( 1C) +Y{ 20 /RE+OZ*Y( 12)-Y( 18) /RE-CZ*Y (9 ) + CX*Y(26) 

Y P ( 22 )=Y (24 )~ORX*Y ( 1 >-TR*Y{ 15 )-OX*Y{ 1C) 

YP (23 ) = Y (25)-CRX*Y (13 )-TR*Y(6 }-OX*Y{ 19 > 

YP ( 24 ) = G*Y( 27 )-CIRX*Y{ 3)-TR*Y( 16)-CX*Y(11 ) 

YP(25 )=-G*Y<26 )-OFX*Y (15 )-TR*Y( 5 )-CX*Y(2C ) 

YP(26)=0RZ*Y (22) + Y<25) /RE+G Z*Y ( 27 > -C RX*Y (17 )-TR*Y { 8 )-CX*Y { 2 1 ) 
YP( 27 )=-Y (24 )/RE-OZ*Y{ 26)+OX*Y{ 28 )-ORX*Y( 1C)-TR*Y( 2C )-CX*Y(12) 
YP<28)=-2.*CPX*=Y<22 ) -2 . * TP*Y ( 25 )-2 . *□ X*Y ( 27 ) 

101 IF (MCDE-2) 1C5» 1C2» IC5 

CGN'T I MGLS FILTER GAINS 

102 C K (.1. ) = Y ( 1 ) / R 
CK(2)=Y(4) /R 
CK ( 3 ) = Y ( 3 ) /R 
CK (4 ).= Y ( 6 )/P 

IF ( A.EG.C. ) GO TC 1C6 
CK(5 )=Y(7)/R 
CK(6)=Y(10)/R 
IF(B.EG.C.) GO TO 106 
CK(7)=Y(13 )/R 
CK ( 8 ) =CK { 7 ) 

CK(S)=Y( 14 )/R 
CK ( 1 0 )=Y ( 15 ) /R 
CK(11) = V( 17) /R 
CK ( 12 )=Y ( 19)/R 
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ck( ib)=yi 2 ? ) /p 

CK { L4)=Y (23)/R 

CONTINUOUS FILTER COMPENSATION TERMS 

ids cTiiisiYmnm+c^Ytnininu/R 
CT(2)={Y(3)*Y{3)+C*Y( H)*Y(Hn/R 
CT ( 3 ) = ( Y ( I )*Y<3 )+C*Y( 1 3 ) *Y { 14) )/R 
CT(4)=(C*Y( 1 3 )*Y ( 1 3 ) 4-Y ( 4 ) *Y (4 ) )/R 
CT(5)=(C rt YI 15 )«v(15) + y(6)*Y(6) )/R 
CT { 6) = IC#Y ( 13)*Y(15)+Y(4)*Y(6) )/R 
IF (A.EQ.O.) GO TO 103 
CT 1 7 ) = ( 8*Y ( 13)*Y(17)+Y{4)*Y(7) )/R 
CT(8)=(B*YI ]5|*Y(17)+Y(6)*Y(7II/R 
CT(9)=(8*Y< 17)*Y(17)+Y( 7)*Y(7) ) /R 
CT( 10)=(Y( 1)*Y( 10)+B*Y( 13)*Y( 19) ) /R 
CT { U> = { YI3)*Y(10) + B*Y{ 14)*Y( 19))'/R 
CT { 12 ) = ( Y ( 10 )*Y( 10)+8*YI 19 > *Y { 19 ) ) /R 
IF IB.EQ.O.) GO TO 103 
CT(13)=(Y(1)*Y( 13)+Y( 13)*Y(4) ) /R 

CT ( 14) = ( Y{ 13 ) #Y( 3) +Y( 4 ) #Y( 14) ) /R 
CT( 15)=( Y( 1 )*Y( 15)+Y{ 13 ) YY{6) ) /R 
CT( 16 ) = ( Y{ 3 )^Y { 15 ) +Y( 14) *Y { 6) ) /R 
CT1 17) = { Y< 1 )*Y{ 17]+YUBI*V(7) ) /R 
CT(18)=(Y(3) *Y (17)+Y{14)*YI7))/R 
CT( 19)=(Y( 13 )*Y(10)+Y(4)*Y( 19) ) /R 
CT { 30 ) - ( Y (1 5 ) *Y ( 1 0 ) +Y I 6 ) *Y ( 19 )) /R 
CTI21 )=(Y( 17)*Y(10)+Y{7)*Y{ 19))/R 
CTI22) ={Y{ 1)*Y( 22)+Y{ 13)*Y<23)) /R 
CT(23)=(Y( 13>*Y{22)+Y<4)*Y(23)>/R 
CT{ ?4)=IY{ 3)*YI 22)+YI 14)*Y(?3) ) /R 
CT I 25 )=(Y{ 15)*Y(22)+Y(6)*Y(23) )/R 
CT( 26)=(Y( 17)*YI22)+Y(7)*Y{23) )/R 
CT ( 27 ) = ( Y( 10)YY( 22J+YI 1 9 ) *Y ( 23 ) ) /R 
CT(28)=(Y(2?)*Y(22)+Y(23)*YI23) )/R 

103 DO 104 1 = 1 , N 

104 YP( l > = YP ( I ) -CT ( I ) 

105 RETURN 
END 


SUBROUTINE KIJTAM 

3(« £ a*#####*### * % ## jjc # # # >[t % # jjnjc % 

INTEGRATION SUBROUTINE; R.KUTT4 4 METHOD, CONTROLLED PRECISION 

DIMENSION YO { 2 8 ) t Yl(28), YP2I28), YP3I28), Y4( 28 ) , Y5( 28) , DERIVI2 
18) 

COMMON A,R,Y(28) ,YP(28) ,R,N, J,NS 

COMMON G »R E , 0X» GZ, ORX , ORZ ,T R , Q , MODE , PAS, E PSI , SEUIL (28) ,T 

LOGICAL S2 

LOGICAL SI 

NS=0 

NIT = 1 
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op 25 1 = 1 , N 
YO(T)=YU) 

25 DFRIV( I )=YP( I) 

S? = . FALSE. 

TO = T 

5 H8 = ° AS/8 . 

H38-P AS#3 . /8. 

HI 5-PAS* 1.8 
H?3=°A$*2./3. 

HO=PAS*2. 

H3= PA S/3. 

H6=PAS/6. 

H 2=PAS / 2 . 

******* ******* ***************** 

INTEGRATION STEPS ; TEST ON POSITIVITY OF MEAN-SQUARED VALUES 


DO 15 1 = 1, N 

Y1 ( n=YO( n+DFR I V { I )*PAS/3 . 
y { n=Yim 
15 CONTINUE 

IF ( Sl( Y( l ) ,Y{4) , Y(2 ) , Y(5 ) ,Y(9) ,YU2) ,Y(?8) } ) GO TO 27 
T = T 0 + H 3 
CALL DAUX 
DO 2 I -I ,N 

2 Y( I)=Y0( I ) + H6U DER I V( I) +YP( I ) ) 

IF(S1 (Yd) ,Y(4J,Y(2 ),YI5 J,Y(9) ,Y( 12) ,Y( 28 m GO TO 27 
CALL DAUX 
DO 4 1=1, N 

4 Y( I ) = YO( d+H8*DERIV{ I ) +H38*Y P ( I ) 

DO 6 1=1, N 
A YP 2 ( I ) =YP ( I ) 

IF<Sl(Ym,Y(4),Y(2),Y(5),Y(9),Y(l?),Y(28))J GO TO 27 

T=T0+H2 

CALL DAUX 

DO 7 1=1, N 

YP3 ( I ) = YP ( I ) 

Y4( I) =YO( I )+H?*DERIV< I )— H15*Y D 2( D+YP3I I ) *HQ 

7 Y(I ) = Y4U) 

IF(S1 IY( 1 ) ,Y(4) ,Y(2 ), Y(5) ,Y (9) ,Y { 12) »Y( 28) ) ) GO TO 27 
T= TO+ P AS 
CALL DAUX 
DO 8 1=1, N 

8 Y5( I) =Y0( I )+H6*DERTV( I ) +YP3 ( I )*H23 + H6*YPI I ) 

IF (Si (Y5(1),Y5( 4),Y5(2)tY5(5),Y5(9) , Y5 ( 1 2) ,Y5(28) ))' GO TO 27 
NTT=N I T+l 

IF (NIT-6) 22,23,23 

23 WRITE (6,24) ER , Y5 ( U , Y 4 { 1 ) 

24 FORMAT U5X,6HP N R ,2X,3E11.4) 

GO TO 19 

****** ************************* 

TEST ON PR FC IS ION 


22 F=0. 

DO 9 1=1, N 

ER=ABS(0.2*(Y5( I >-Y4( I ) ) /AMAX1 ( ABS (Y5( I ) ) ,SEUIL( I) ) ) 
IF (ER.GT.E) E= ER 
9 CONTINUE 
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IF (E-EPSI) 13,13,16 

* * # * # # # afr # # *e £ # * * Jfe ^ ^ # # # * # >*t # * # A £ ={c 

IF PRECISION NOT RE ACHED, ST &> SIZE HALVED 


16 P AS=H2 
GO Tn 5 

13 IF (64.*E— EPSI I 18,19,19 
19 $2=.TRUE. 

18 DO 29 1 = 1, N 

29 Y ( I ) = Y5I I) 

CALL DAUX 
GO TO 26 

27 NS=NS+1 

WRITE (6,30) NS,T, Y( 1 ),Y (2) ,Y(4) ,YI 5) ,Y(9) ,Y( 12) ,Y(28) 

30 FORMAT (I4,2X,8E11.4) 

# jfc i}s £#5)c5)s>(£:$c a{c^s2{<3}c^eajs5{«a(s55e3(e^e^c^e 

IF POSITIVITY TEST CALLED TOO MANY TIMES,STOP 

IF (NS.GE.5Q) GO TO 21 
PAS=P AS/2. 

T=TO 
GO TO 5 

************************* ****** - 

IF PRECISION EXCESSIVE, STEP SIZE DOUBLED 

! 26 IF ( S2 1 GO TO 21 

P AS=HD 
21 RETURN 
END 


SUBROUTINE UPDAT 
* * * * # * * * * * * * * * * * * * * * * * * * * * * * ** * 

TO COMPUTE DISCRETE FILTER GAINS AND UPDATE COVARIANCE MATRIX 

DIMENSION HK ( 1 4 ) 

COMMON A,8»Y{28) ,YP(28) ,R,N,J,NS 
DO 200 1=1,14 

200 HK( I ) =0. 

J = 1 

C = A*B 

AK 1=Y ( 1) +R 
AK2=Y ( 4)+R 

0=AK1*AK2-C*Y( 13)*Y(13) 

IF CD) 201,201,202 

201 PRINT 1201 

1201 FORMAT { 25X, 15H INF IN ITE GAINS) 

J=J+i 
GO TO 204 

* * * * * * * * * * * * * * * * # * * * * * * * * * * * * * * 

DISCRETE FILTER GAINS 

202 HK(l)r CY I 1 ) *AK2-C*Y { 1 3 ) *Y ( 13 )) /D 
HK(2)=<-C*Y(13)*Y{ 13) +Y I4)*AK1)/D 
HK < 3) = { Y { 3 ) *AK2-C*Y ( 1.3-) *Y ( H) ) /D 
HK ( 4) = {— C#Y ( 13 ) *Y ( 15) +Y(6)*AK1 )/D 
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IF (A. EQ.O. ) GO TO 203 
HK ( 5) = C-B*Y (13)*YI 1 7 ) +Y { 7 } *AK1 ) /D 
HK(6)=(Y{ 10)*AK2-R*YI 13I*Y( 19) >/D 
IF (3. EQ.O.) GO TO 203 
HK< 7 ) = {-YU )*Y(13)+Y< 13)*AK1 )/D 
HK t Q) = C Y { 13)*AK2-Y( 13)*Y(4))/0 
HK ( 9 ) ■=■( — Y ( 1 3 ) * Y ( 3 ) +Y ( 14)*AK1)/D 
HK ( 10 ) = C Y( 15)*AK2-Y(13)*Y{6) )/D 
HK { 11 )=( Y{ 17)*AK2-Y( 13)*Y< 7) ) /D 
HK( 12 ) = { — Y < 13)*Y(10)+Y(19)*AK1 )/0 
HK (13 ) - l Y( 22 )#AK2— Y { 13) #YI 23 ) ) /D 
HK( 14 )=C-Y{ 13)*Y(22)+Y(23)*AK1)/D 
$*#*$* **#*## :**£*;$:£ ******* 

UPDATE COVARIANCE MATRIX 

203 Y( 1 ) = Y ( 1 ) -HK (1 ) *Y(1)-C*HK< 7 )*Y( 13) 

Y( 2)=Y( 2)-HK{3)*Y{3 )-C*HK(9 )*Y( 14) 
Y(3)=Y(3) -HK( 1)*Y( 3)-C*HK{ 7)*Y{ 14) 
Y(4)=Y(4)-HK(2 ) *Y ( 4 ) -C*HK { 8 ) *Y { 13) 
Y(5)-Y(5)-HK(4) *Y ( 6 )-C*HK C 10 ) *Y ( 15) 

Yf 6)=Y(6)-HK(2)*Y{6)-C*HK<8)*Yf 15) 

IF { A . EQ.O . ) GO TO 204 

Y( 7)=Y(7)-HK(2)*YI7)-B*HK(8)*Y( 17) 

Y ( 8 ) = Y ( 8 )-HK(4)#Y(7 ) -B*HK { 1 0 ) *Y f 17) 

Y( 9)=Y(9)-HK{5) *Y(7)-B*HK( 11 )*Y( 17) 

Y{ 10) =Y( 10 ) -HK { 1 )*Y( 10 )— B*HK ( 7) *Y { 19) 
Y{ ll)=YUl)-HK{ 31*Y( 10)-B*HK(9)*Y< 19* 

Y ( 12)-Y(12)-HK(6)*YI10 ) -8*HK ( 12 ) *Y( 19) 

IF (8. EQ.O.) GO TO 204 

Y( ] 3) =YI 13)-HK( l)*Y(13)-HKf 7)*Y{4) 

Y I 14) =Y { 14 ) - HK { 2 ) *Y (14) -HK ( 8 ) *Y t 3) 

Y ( 15)=Y( 15.)-HK (1)*YU5)-HK{7)*Y(6) 

Y ( 16) =Y{ 16)— HK(3)*Y(15 )— HK ( 9) *Y( 6 ) 
Y(17)=Y( 17)-HK(1)*Y< 17)-HK{ 7}*Y{?) 

Y( 18) = Y{ 18 > -HK (3)*YI17)-HK{9)*Y{7) 

Y { 19 ) =Y ( 19 }— HK( 2 )*Y( 19 )-HK f 8 ) *Y ( 10 ) 
Y(20)=Y{20)-HK(4)*YI19 l-HKUO)n(IO) 

Y( 21) =Y<21)-HK{5)*Y(19)-HKI11)*YC 10) 
Y(22)=Y(22)-HK{ 1 ) *Y ( 2 2 )-HK< 7 ) *Y{ 23) 

Y( 23) = Y{ 23 )— HK (2 ) *Y (23 ) — HK ( 8 ) #Y( 22) 
Y(24)=Y(24)-HK(3)*Y(22)-HK(9)*Y(23) 

Y { 2 5) =Y ( 2 5 ) -HK ( 4 ) #Y { 23 ) — HK ( 10)*Y(22) 

Y I 26)=YI26)-HKI 5)*Y<23 ) -HK C 1 1 ) *Y { 2 2 ) 

Y( 27)=Y{27)-HK(6)*Y< 22)-HK( 12)*Y(23) 
Y{28)=Y(28)-HK(13)*Y(2?)-HK<14)*Y(23) 

204 CALL DAUX 
RETURN 
END 


LOGICAL FUNCTION S 1 1 A , B, C, 0* E*F, G ) 

■’ll * * X-- 4 =* # # # # ## £ # £ * * £ $ A $ 3fc Jft ifr ifc :$( * # $ 

TO CHECK IF THE DIAGONAL TERMS OF CDV. MATRIX ARE POSITIVE 

IF (G.LT.O . .AND.G.GT.-l. E— 1 B » G=+0.0 
Sl=. FALSE. 

IF{ A.LT.O..QR.B.LT.O. .OR.C.LT.O. .OR. D.LT .0. .OR. E .LT. 0. .OR.F.LT.O.. 
10R. G.LT.O.) Sl= .TRUE. 

RETURN 

END 
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TO PLOT RMX AND RVX IN MODES 1,2 AND 3 WITH MODELS 1 ,2 
AND 3 DIFFERENT OPT 


DIMENSION TI (50) ,XU( 50>,X21(50) ,X31(50> ,X12(50) ,X?2(50) ,X32(50),V 
lll<50> ,V?1( 50), V31(50), V12I50I ,V?2< 50) , V32 ( 50 ) , XI 31 ( 50 ) , X132(50 ) ,X 
2133(50) f X231( 50) , X232 ( 70 ) ♦ X233 ( 50) , V131 ( 50 ) , VI 32 ( 50 ) ,V133 ( 50 ) T V231 
3(50)»V232( 50),V233(50) 

COMMON AQ, AR ,OPT1,OPT2,OPT3 

CALL NEWPLT( » M6 175 •,» 7312* * 1 WHITE ' * "BLACK* ) 
******************************* 

READ THE MAX. NUMBER OF CASES TO BE PLOTTED 

1010 FORMAT ( 13 ) 

READ (5,1010) K2 
******************************* 

READ RESULTS OF THE PREVIOUS PROGRAM ON CARDS 
CHECK THE ORDER OF THE CARDS 

1 READ (5,1020) K,AQ,AR 

1020 FORMAT (5HCASE I2,2X,4HN = Elt.4,3X,4HR = Ell. 4) 

K 1 =K 

READ (5,1021) OPTl,OPT2,CPT3 

1021 FORMAT (7H0PT1 = E l 1. 4, 3X , 7H0PT2 = Ell .4, 3X,?H0PT3 = Ell. 4) 

0PT1=0PT 1/60 . 

0PT2=0PT2/60. 

0PT3=0PT3/60. 

READ (5,1022) X M 1 , XM2 , VM l ,V M2 

1022 FORMAT ( 6HXM1 = E11.4,6HXM2 = E11.4,6HVM1 ® E11.4,6HVM2 = Ell. 4) 
READ (5,1023) XM 13, XM23, VM13, VM23 

1023 FORMAT 17HXM13 = E11.4,7HXM23 = Elt.4,7HVM13 = £li.4,7HVM23 = Ell. 4 
l) 

DO 100 M=l,50 

READ (5,1024) T I ( M ) , Xll ( M) , X21 ( M ) , X3 1 ( M ) ,X 1 2 ( M) , X22 ( M ) , X32 ( M ) , VI 1( 
1M),V21(M),V31(M) ,V12(M),V22(M),V32(M) ,K,M1 
IF (K.NE.K1.0R.M1.NE.M) GO TO 1 

READ (5,1024) T I( M ) , X13U M I ,X 132 (M ) ,X 1 33( M ) , X 23 l( M ) , X232 ( H ) ,X233(M 
1),V131 (M) ,V132(M) ,V133(M) ,V231(M) ,V232(M) , V233( M ) , K , Ml 
IF (K.NE.K1.0R.M1.NE.M) GO TO 1 

1024 FORMAT (F6.0, 12F5.3,6X*I2*2X,I3) 

100 CONTINUE 

******************************* 

COMPUTE REAL VALUES WITH THE NORMALIZATION CONSTANTS 

DO 101 M=1,5C 

IF (X2 1( M) . L-T , XM23) NPl=M+5 
IF(V21(M).LT.VM23) NP2=M*5 
TI ( Ml =T H M) /60 • 

XU(M) = X11(M)*XMI 
X21(M)=X21(M)*XM1 
X3l(M)=X31(M)*XMl 
X12(M)=X12(MI*XM2 
X22(M)=X22(M)*XM2 
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X32{M)=X32(M)*XM2 
Vll(M|=Vll { M)*VM1 
V2HM )=V21( M)*VM1 
V31{M)=V31(*)*VM1 
V12(H)=V12(H)*VH2 
V22(M )=V22( MKV«2 
„ V32(M)=V32( M)*VM2 

X131( M) = X131 (HI *XM 1 3 
* X132( H)=X132( M)*XM13 

X 1 33( M)=X133(M)*XH13 
X231( M )=X231(M)*XM23 
' X232( M)=X2?2(M) *XH23 

X233{M)=x 233(H) *XM23 
V13K M)=V131 (H) *VM13 
V 132( M ) =V1 3 2( M I *VMi3 
V 1 33 ( H I =V13 3 ( H ) * VM 13 
V231(M)=V 231(H) *VH23 
V232(MI=V232(M)*VM23 

101 V233( M)=V233(«)*VM23 
**##$*****#******************** 

PLOT ALL THESE RESULTS 

CALL IOENPL (1,0) 

CALL PICTUR ( 6.»4.,*T HI N» , 5 , »RMX • ,3,TI f X 11 , 50 ,0 . » KS , TT , X2 1, 50 ,0 . , K 
1S,TI, X31,50,0.,KS) 

CALL IOENPL (1,0) 

CALL PICTUR <6., 4. ,*T MIN* , 5, • RVX • , 3 ,TI ,V 11 , 50 ,0. , KS, TI , V21 , 50,0. , K 
1S,TI,V31,50,0.,KS) 

CALL IOENPL (3,1) 

CALL PICTUR (6., 4. , *T MIN* ♦ 5, «RHX* ,3,TI , X12,50 ,0. , KS,TI , X131 ,50,0., 
1KS,TI ,X 132, 50,0. , KS ,T I , X 133, 50 ,0 . , KS) 

CALL I DENPL ( 3*1.) 

CALL PICTUR(6.,4.,*T HIN ■ , 5 » *RVX • , 3, T I ,V 12 » 50 , 0. , K S, TI , VI 31 , 50 ,0 . , 
1KS,TI ,V132,50,0.,KS,TI,V133,50,0.,KS) 

CALL IOENPL (3,2) 

CALL PICTUR (6.,4.,*T MIN*, 5, •RHX* ,3,TI ,X 22 , 50 , 0 . , KS, TI , X231 , 50, 0. , 
IK S, TI ,X232,50,0. , KS ,T I ,X233 , 50 , 0. , KS) 

CALL I0ENPL(3,2 ) 

CALL PICTUR (6, ,4. ,*T HI N» • 5, • RVX* , 3, TI , V2? , 50 ,0 . , KS » TI , V23 I , 50, 0 . » 
IKS ,T I ,V232,50,0. f KS,TI,V233,50,0.,KSI 
DO 10 2 H-NP 1,50 

102 X21LM) = 0. 

DO 103 M=NP2, 50 

103 V2 l ( H | =0 . 

CALL IDE NPL( 0,2) 

CALL PICTUR (6.*4., , T M IN* « 5, • RMX» ,3 ,TI , X 21 » 50 ,0 . t KS ,TI , X22 , 50 ,0. ,K 
1S,TI,X233,50,0.,KS,TI,X231,50,0.,KS) 

CALL IDENPL(0,2) 

CALL PICTUR (6.,4.»*T H IN • , 5, • RVX* , 3, TI , V21 , 50 ,0. , KS ,T I , V22 . 50,0. , K 
IS,TI,V233,50,0. ,KS,TI, V231,50,0. ,KS) 

* **************** ************** 

IF NUMB. OF CASES PLOTTED = NUMB. DESIRED , STOP 

IF (K1.EQ.K2) GO TO 2 
GO TO 1 
2 CALL ENDPLT 
' STOP 
END 
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SUBROUTINE IDENPL(J,L> 

#«**«***#*(***>!<*$ *#***« *«***$**$ 

TO DRAW SOME IDENTIFICATION ON TOP OF EACH PLOT 

COMMON AQ» AR ,OPT 1 ,OPT 2 » 0PT3 
CALL PLOTU .5,3.25,-31 
CALL PLOTUO. ,0.75,-2) 

CALL PLOTU 1.5,0., -2) 

CALL PLOTUO. ,-0.75, -2) 

CALL PL0TlI-1.5,0.,-2) 

IF(J.NE.O) GO TO 2 

CALL SYMBL5(0.1»O.5,.2,' MODEL* «,0.,47) 

CALL NUMBRl(1.3,0.5,.2,Lt0.,-l) 

GO TO 3 

2 CALL SYM3L5(0.1,0.5,.2,*MODE * »,0.,*7) 

CALL NUMBR1(1.3,0.5,.2,J,0.,-1I 

3 CALL SYMBL5 (0.05,0.2*.05»*N = »,0.,+4) 

CALL NUMBR1 (O.2»0.2». 10, AQ,0«, fl-2 ) 

CALL SYM8L5 (0.05, .05, .05, *R « »,0.,+4) 

CALL NUM3R1(0.2,.05,. 10,A*,0.,*0) 

IF(J.EQ.O) GO TO 4 

CALL SYMBL5 (0.9 ,0.35, .05, •0PT1 ■= »,0.,+7) 

CALL NUMSR 1(1.2, 0.35, .10.0PTl,0 .,+-1) 

CALL SYMBL5(0.9,0.2,.05, * 0PT2 = *,0.,+ 7) 

CALL NUMBRl (1,2 ,0.2 , .10, 0PT2,0 . ,+1) 

4 CALL SY M BL 5(0. 9,0.05,. 05**0 PT 3 = *,0.,+7) 

CALL NUMBRl ( 1.2, 0.05, .10, DPT‘3,0. ,+ l> 

IF (L.EQ.O) GO TO 1 
IF(J.EQ.O) GO TO 5 

CALL SYM9L5 (.05,. 35, .05, • MODEL = »,0.,+8) 

CALL NUMBRl (0.50, .35, . 10*L,0.,-1 I 
GO TO 1 

5 CALL SVMBL5 (0.05,. 35,. 10, ‘MODE l 2 3*,0.,+10) 
l CALL PLOTU— 0.5»— 3.25,— 3) 

RETURN 

END 
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